{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 18\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from the previous chapter\n",
    "\n",
    "Read the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('data/glucose_insulin.csv', index_col='time');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolate the insulin data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<function modsim.modsim.interpolate.<locals>.wrapper(x)>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = interpolate(data.insulin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1e-05"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G0 = 290\n",
    "k1 = 0.03\n",
    "k2 = 0.02\n",
    "k3 = 1e-05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate basal levels, we'll use the concentrations at `t=0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gb = data.glucose[0]\n",
    "Ib = data.insulin[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the initial condtions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>G</th>\n",
       "      <td>290</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>X</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "G    290\n",
       "X      0\n",
       "dtype: int64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init = State(G=G0, X=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "182"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_0 = get_first_label(data)\n",
    "t_end = get_last_label(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>G0</th>\n",
       "      <td>290</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k1</th>\n",
       "      <td>0.03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k2</th>\n",
       "      <td>0.02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k3</th>\n",
       "      <td>1e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>G    290\n",
       "X      0\n",
       "dtype: int64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gb</th>\n",
       "      <td>92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ib</th>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>I</th>\n",
       "      <td>&lt;function interpolate.&lt;locals&gt;.wrapper at 0x7f...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "G0                                                     290\n",
       "k1                                                    0.03\n",
       "k2                                                    0.02\n",
       "k3                                                   1e-05\n",
       "init                        G    290\n",
       "X      0\n",
       "dtype: int64\n",
       "Gb                                                      92\n",
       "Ib                                                      11\n",
       "I        <function interpolate.<locals>.wrapper at 0x7f...\n",
       "t_0                                                      0\n",
       "t_end                                                  182\n",
       "dt                                                       2\n",
       "dtype: object"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = System(G0=G0, k1=k1, k2=k2, k3=k3,\n",
    "                init=init, Gb=Gb, Ib=Ib, I=I,\n",
    "                t_0=t_0, t_end=t_end, dt=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Updates the glucose minimal model.\n",
    "    \n",
    "    state: State object\n",
    "    t: time in min\n",
    "    system: System object\n",
    "    \n",
    "    returns: State object\n",
    "    \"\"\"\n",
    "    G, X = state\n",
    "    k1, k2, k3 = system.k1, system.k2, system.k3 \n",
    "    I, Ib, Gb = system.I, system.Ib, system.Gb\n",
    "    dt = system.dt\n",
    "    \n",
    "    dGdt = -k1 * (G - Gb) - X*G\n",
    "    dXdt = k3 * (I(t) - Ib) - k2 * X\n",
    "    \n",
    "    G += dGdt * dt\n",
    "    X += dXdt * dt\n",
    "\n",
    "    return State(G=G, X=X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "        \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    t_0, t_end, dt = system.t_0, system.t_end, system.dt\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_simulation(system, update_func);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical solution\n",
    "\n",
    "In the previous chapter, we approximated the differential equations with difference equations, and solved them using `run_simulation`.\n",
    "\n",
    "In this chapter, we solve the differential equation numerically using `run_ode_solver`, which is a wrapper for the SciPy ODE solver.\n",
    "\n",
    "Instead of an update function, we provide a slope function that evaluates the right-hand side of the differential equations.  We don't have to do the update part; the solver does it for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes derivatives of the glucose minimal model.\n",
    "    \n",
    "    state: State object\n",
    "    t: time in min\n",
    "    system: System object\n",
    "    \n",
    "    returns: derivatives of G and X\n",
    "    \"\"\"\n",
    "    G, X = state\n",
    "    k1, k2, k3 = system.k1, system.k2, system.k3 \n",
    "    I, Ib, Gb = system.I, system.Ib, system.Gb\n",
    "    \n",
    "    dGdt = -k1 * (G - Gb) - X*G\n",
    "    dXdt = k3 * (I(t) - Ib) - k2 * X\n",
    "    \n",
    "    return dGdt, dXdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-5.9399999999999995, 0.0)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope_func(init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run the ODE solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2, details = run_ode_solver(system, slope_func, t_eval=data.index);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`details` is a `ModSimSeries` object with information about how the solver worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`results` is a `TimeFrame` with one row for each time step and one column for each state variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>G</th>\n",
       "      <th>X</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>290</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>278.476</td>\n",
       "      <td>0.00015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>267.464</td>\n",
       "      <td>0.00147812</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>255.899</td>\n",
       "      <td>0.00330258</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>244.424</td>\n",
       "      <td>0.00428352</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>233.421</td>\n",
       "      <td>0.0048796</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>222.905</td>\n",
       "      <td>0.00539312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>212.909</td>\n",
       "      <td>0.00580811</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>203.454</td>\n",
       "      <td>0.00610843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>194.566</td>\n",
       "      <td>0.00630605</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>186.24</td>\n",
       "      <td>0.00643892</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>178.459</td>\n",
       "      <td>0.00655891</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>171.185</td>\n",
       "      <td>0.0066622</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>164.396</td>\n",
       "      <td>0.00673793</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>158.068</td>\n",
       "      <td>0.00679316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>152.172</td>\n",
       "      <td>0.00686423</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>146.67</td>\n",
       "      <td>0.00695603</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>141.532</td>\n",
       "      <td>0.00703975</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>136.742</td>\n",
       "      <td>0.00708884</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>132.286</td>\n",
       "      <td>0.00710463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>128.148</td>\n",
       "      <td>0.00708845</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>124.315</td>\n",
       "      <td>0.00704154</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>120.769</td>\n",
       "      <td>0.00696712</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>117.496</td>\n",
       "      <td>0.00686816</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>114.481</td>\n",
       "      <td>0.00674565</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>111.709</td>\n",
       "      <td>0.0066005</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>109.167</td>\n",
       "      <td>0.0064336</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>106.84</td>\n",
       "      <td>0.00625981</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>56</th>\n",
       "      <td>104.71</td>\n",
       "      <td>0.00609282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>58</th>\n",
       "      <td>102.762</td>\n",
       "      <td>0.00593238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>124</th>\n",
       "      <td>86.9662</td>\n",
       "      <td>0.00107406</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>126</th>\n",
       "      <td>87.0883</td>\n",
       "      <td>0.000956516</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>128</th>\n",
       "      <td>87.2224</td>\n",
       "      <td>0.000845541</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>130</th>\n",
       "      <td>87.3667</td>\n",
       "      <td>0.000740875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>132</th>\n",
       "      <td>87.5196</td>\n",
       "      <td>0.000642273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>134</th>\n",
       "      <td>87.6796</td>\n",
       "      <td>0.000549496</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>136</th>\n",
       "      <td>87.8454</td>\n",
       "      <td>0.000462316</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>138</th>\n",
       "      <td>88.0157</td>\n",
       "      <td>0.000380513</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>140</th>\n",
       "      <td>88.1896</td>\n",
       "      <td>0.000303877</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>142</th>\n",
       "      <td>88.3658</td>\n",
       "      <td>0.000232205</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>144</th>\n",
       "      <td>88.5436</td>\n",
       "      <td>0.000164302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>146</th>\n",
       "      <td>88.7224</td>\n",
       "      <td>9.90618e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>148</th>\n",
       "      <td>88.9018</td>\n",
       "      <td>3.63786e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>89.0813</td>\n",
       "      <td>-2.38474e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>152</th>\n",
       "      <td>89.2606</td>\n",
       "      <td>-8.17126e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>154</th>\n",
       "      <td>89.4392</td>\n",
       "      <td>-0.000137309</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>156</th>\n",
       "      <td>89.617</td>\n",
       "      <td>-0.000190727</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>158</th>\n",
       "      <td>89.7936</td>\n",
       "      <td>-0.00024205</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>160</th>\n",
       "      <td>89.9687</td>\n",
       "      <td>-0.000291362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>162</th>\n",
       "      <td>90.1422</td>\n",
       "      <td>-0.000338741</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>164</th>\n",
       "      <td>90.3138</td>\n",
       "      <td>-0.000385262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>166</th>\n",
       "      <td>90.4837</td>\n",
       "      <td>-0.00043192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>168</th>\n",
       "      <td>90.6521</td>\n",
       "      <td>-0.000478709</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>170</th>\n",
       "      <td>90.8191</td>\n",
       "      <td>-0.000525623</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>172</th>\n",
       "      <td>90.9848</td>\n",
       "      <td>-0.000572659</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>174</th>\n",
       "      <td>91.1493</td>\n",
       "      <td>-0.00061981</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>176</th>\n",
       "      <td>91.3128</td>\n",
       "      <td>-0.000667074</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>178</th>\n",
       "      <td>91.4754</td>\n",
       "      <td>-0.000714445</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>180</th>\n",
       "      <td>91.6372</td>\n",
       "      <td>-0.000761918</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>182</th>\n",
       "      <td>91.7983</td>\n",
       "      <td>-0.000809491</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>92 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           G            X\n",
       "0        290            0\n",
       "2    278.476      0.00015\n",
       "4    267.464   0.00147812\n",
       "6    255.899   0.00330258\n",
       "8    244.424   0.00428352\n",
       "10   233.421    0.0048796\n",
       "12   222.905   0.00539312\n",
       "14   212.909   0.00580811\n",
       "16   203.454   0.00610843\n",
       "18   194.566   0.00630605\n",
       "20    186.24   0.00643892\n",
       "22   178.459   0.00655891\n",
       "24   171.185    0.0066622\n",
       "26   164.396   0.00673793\n",
       "28   158.068   0.00679316\n",
       "30   152.172   0.00686423\n",
       "32    146.67   0.00695603\n",
       "34   141.532   0.00703975\n",
       "36   136.742   0.00708884\n",
       "38   132.286   0.00710463\n",
       "40   128.148   0.00708845\n",
       "42   124.315   0.00704154\n",
       "44   120.769   0.00696712\n",
       "46   117.496   0.00686816\n",
       "48   114.481   0.00674565\n",
       "50   111.709    0.0066005\n",
       "52   109.167    0.0064336\n",
       "54    106.84   0.00625981\n",
       "56    104.71   0.00609282\n",
       "58   102.762   0.00593238\n",
       "..       ...          ...\n",
       "124  86.9662   0.00107406\n",
       "126  87.0883  0.000956516\n",
       "128  87.2224  0.000845541\n",
       "130  87.3667  0.000740875\n",
       "132  87.5196  0.000642273\n",
       "134  87.6796  0.000549496\n",
       "136  87.8454  0.000462316\n",
       "138  88.0157  0.000380513\n",
       "140  88.1896  0.000303877\n",
       "142  88.3658  0.000232205\n",
       "144  88.5436  0.000164302\n",
       "146  88.7224  9.90618e-05\n",
       "148  88.9018  3.63786e-05\n",
       "150  89.0813 -2.38474e-05\n",
       "152  89.2606 -8.17126e-05\n",
       "154  89.4392 -0.000137309\n",
       "156   89.617 -0.000190727\n",
       "158  89.7936  -0.00024205\n",
       "160  89.9687 -0.000291362\n",
       "162  90.1422 -0.000338741\n",
       "164  90.3138 -0.000385262\n",
       "166  90.4837  -0.00043192\n",
       "168  90.6521 -0.000478709\n",
       "170  90.8191 -0.000525623\n",
       "172  90.9848 -0.000572659\n",
       "174  91.1493  -0.00061981\n",
       "176  91.3128 -0.000667074\n",
       "178  91.4754 -0.000714445\n",
       "180  91.6372 -0.000761918\n",
       "182  91.7983 -0.000809491\n",
       "\n",
       "[92 rows x 2 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the results from `run_simulation` and `run_ode_solver`, we can see that they are not very different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f0d543fc240>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD9CAYAAAC/fMwDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3zU9Z3v8ddvrpmEW4BwEQIo0i9gkQqoq7aSirrFWlu3xW6rRQ/20cupPWfrto9dT7u6tiprt+2pSh+Htmsv2LqslFpFxepaoaLWclcBv3gBDfeEXEnmPnP++E1wCEkIJOSXzLyfj0eYzO8y+cwv4T3f+f6+8/052WwWEREpfD6vCxARkb6hwBcRKRIKfBGRIqHAFxEpEgp8EZEiEfC6gM4YY8LA+cB+IO1xOSIiA4UfGAust9bG81f028DHDfsXvC5CRGSA+giwLn9Bfw78/QC//e1vGTNmjNe1iIgMCAcOHOD666+HXIbm68+BnwYYM2YM48eP97oWEZGB5riucJ20FREpEgp8EZEiocAXESkS3erDN8ZcDdwDnAkcAr5vrf1pbuhkM5DI2/wla+2Vuf2uy+03FlgL3GStPdSL9YuISDedMPCNMWOB3wHXWmtXG2NmAS8aY9bjvkOos9YeN4zGGDMdeBCYD2wA7gWWA5f1Yv1dWrOxmmWrd1BbH2VkeYSF86dRNbuyr368iEi/csIuHWvtfqAiF/Y+YASQwm3Zzwa2dLLrDcAqa+06a20MuA24xBgzpXdK79qajdUsWbGVmvooWaCmPsqSFVtZs7G6L368iEi/060+fGttszGmFIgDzwA/sda+CcwCRhljXjXGHDTGrDDGjMvtNh3YnvcYrUA1MKNXn0Enlq3eQTx57KikeDLNstU7+uLHi4j0Oydz0jYGlOF+AnaRMeZmoAV4EZgHGCAKPJrbfhDQ2u4xWoHSnhTcXbX10ZNaLiJS6Lr9wStrbQb35OwGY8zPgE9aa6/J38YYcytQY4ypxH0xiLR7mFLgSM9K7p6R5RFqOgj3keXtSxIRKQ4nbOEbY+YaYza2WxwGGowx3zXGTMtbHsrdxnC7c0ze45QCE8jr5jmdFs6fRjjoP2ZZOOhn4fxpnewhIlLYutPC3wKMy7Xe7wMuBG4GrgW+Acwxxnw+t+19wJPW2hpjzMPAOmNMFfAysBjYbK3d2cvPoUNto3E0SkdExHXCwLfWNhpjrgLuB+7APfH6RWvtWmPM67nlb+Ue60ngS7n9XjPGLAKWAuOAV4AFp+VZdKJqdqUCXkQkp1t9+NbaTcCHO1h+GLi+i/1WAitPuToREek1mlpBRKRIKPBFRIqEAl9EpEgo8EVEioQCX0SkSCjwRUSKhAJfRKRIKPBFRIqEAl9EpEgo8EVEioQCX0SkSCjwRUSKhAJfRKRIKPBFRIqEAl9EpEgo8EVEioQCX0SkSCjwRUSKhAJfRKRIKPBFRIqEAl9EpEgo8EVEioQCX0SkSCjwRUSKhAJfRKRIKPBFRIpEoDsbGWOuBu4BzgQOAd+31v7UGBMClgCfAdLAj6y1i/P2uy6331hgLXCTtfZQ7z4FERHpjhO28I0xY4HfAf9krR0MLAB+bIyZBdwJGGAycD5wozFmYW6/6cCDwE3ACOBNYPlpeA4iItINJwx8a+1+oMJau9oY48MN7xTQDNwI3G2trbfW7gZ+AHw5t+sNwCpr7TprbQy4DbjEGDPlNDwPERE5gW714Vtrm40xpUAceAb4CVCD21WzPW/TN4AZue+n56+z1rYC1XnrRUSkD3WrDz8nBpQB5wJPAdHc8ta8bVqB0tz3g9qta79eRET6ULcD31qbARLABmPMz4A5uVWRvM1KgSO571varWu/XkRE+lB3TtrONcZsbLc4DNQDB3BP2raZyvvdONvz1+W6hCZwbBeQiIj0ke608LcA44wxtwL3ARcCNwPX4gb+HcaYV3G7cL6Z2wbgYWCdMaYKeBlYDGy21u7s1WcgIiLd0p1ROo3AVcDfAXXAz4AvWmvXArcDrwPbgPXASmBpbr/XgEW5+7XAObhDOkVExAPd6sO31m4CPtzB8hjwtdxXR/utxH0REBERj53MKJ0BY+0vlpCNt3Dpl7+Fz/f+m5g1G6tZtnoHtfVRRpZHWDh/GlWzKz2sVESk7xTkXDql+zcxoe6v2L/+5eiyNRurWbJiKzX1UbJATX2UJSu2smZjtXeFioj0oYIM/JbRHwLg0CtPH122bPUO4sn0MdvFk2mWrd7Rp7WJiHilIAN/yrxryGRhbPMOmusOA1BbH+1w286Wi4gUmoIM/HFnTmJvcBIBJ8O2Z1cBMLK8/WfA6HK5iEihKcjAByiZMQ+A4NsvkM1mWTh/GuGg/5htwkE/C+dP86I8EZE+V7CBP3Pe5TRmSinPNrBr01+pml3JLQtmUlEewQEqyiPcsmCmRumISNEoyGGZACXhELWj5jC09s/sf+kpzpp9IVWzKxXwIlK0CraFDzC56moyWaho3E60sd7rckREPFXQgX+Wmcy7/okEnAzbn33c63JERDxV0IEPEPzg5QD4d64hm0mfYGsRkcJV8IE/5/LLOZwZxJBsE7vWr/O6HBERzxR84JdGQhwecxEANS8/6XE1IiLeKfjAB5h25TUksn4qWt6maf97XpcjIuKJogj8iRPPYFd4KgD2j7/3uBoREW8UReADDL/wKgDK9v6VdFzz54hI8SmawJ99yQVUZ0ZRQpyda58+8Q4iIgWmaAI/4PcRPetSAFq3PEM2m/W4IhGRvlU0gQ8w+2Mf50gmTHnyEAd3bve6HBGRPlVUgV8xYgj7hs4E4O3n9clbESkuRRX4AJOqPgHA8NottDY3eVyNiEjfKbrANzOmsdc3jrCTYusf9UEsESkeRRf4juMQ/uBl7p2da3XyVkSKRtEFPsCHrvgYrdkwo7M1vP7XjV6XIyLSJ4oy8EMlJTSNmQPAvhfVrSMixaEoAx/AXPEpACa0bmf//lqPqxEROf26dYlDY8wVwL8BU4BDwL9ba39qjAkDzUAib/OXrLVX5va7DrgHGAusBW6y1h7qxfpP2Wu1Qf6jcQGN6RIG3/dnvrTgfF3+UEQK2gkD3xhTCawEbgQeA2YDfzTG7AYOA3XW2jEd7DcdeBCYD2wA7gWWA5f1Uu2nbM3Gapas2Eo8HQGgOelnyYqtAAp9ESlY3enSmQQ8bK191FqbsdauB9YAl+CG/5ZO9rsBWGWtXWetjQG3AZcYY6b0vOyeWbZ6B/HksVe/iifTLFu9w6OKREROvxO28K21LwAvtN03xgwHPgI8BHwMGGWMeRUYDfwZ+Adr7V5gOm7Lvu1xWo0x1cAM4M3efBInq7a+49kyO1suIlIITuqkrTFmKPA48Apu904L8CIwDzBAFHg0t/kgoLXdQ7QCpT2ot1eMLI90uHzY4FAfVyIi0ne6ddIWwBjzAdyQ3w5cb63NALe22+ZWoCbX798CtE/WUuBIjyruBQvnT3P78PO6dYKkOG+QploQkcLVrRa+MeZS3Fb9H4DP5PrkMcZ81xgzLW/TtiZyDPeFweQ9RikwIbfcU1WzK7llwUwqyiM4QHnE4e/LXuKC6FpaookT7i8iMhB1Z5TOZOAJ4NvW2gfarT4XmGOM+Xzu/n3Ak9baGmPMw8A6Y0wV8DKwGNhsrd3Za9X3QNXsyqMjcrLpJNvvfYRINsorz7/IZVd91OPqRER6X3e6dL4GDAYWG2MW5y3/CXAzcD/wVu6xngS+BGCtfc0YswhYCozDfYewoPdK7z1rtxzgP5qvpTHhZ9DzNTijqvnoHA3PFJHC0p1ROrfSrq++neu72Hcl7hj+fuvomPykeyiOZMI88MhmHEdj8kWksBTt1AptOhqTn0xnNSZfRApO0Qd+Z2PvazQmX0QKTNEHfmdj8geHNE++iBSWog/8hfOnEQ76j1kWJMXFgW2k0xmPqhIR6X1FH/jtx+RXDC3hM2WvcFXJJjb9davX5YmI9Jpuf9K2kOWPyQdY/+AWOPA2B196Ai46z8PKRER6T9G38Dty9hWfBuCs6Ou8V33Q42pERHqHAr8D5RPOpjYyiRInxetPP3riHUREBgAFfidGffiT7u2Bl2k+oiGaIjLwKfA7MfH8D9PoG8Zw3xHWP/201+WIiPSYAr8TjuPD/8ErAPDb/yad0bh8ERnYFPhdmHbFNcQIUckBNq17xetyRER6RIHfhWBJKc3j/gaA+ldWeVyNiEjPKPBPYNr860hnHc6MW97ZucvrckRETpkC/wSGjB5L7dDp+J0sb/7xEa/LERE5ZQr8bph4xWcBqGzcRM2BWo+rERE5NZpaoQNrNlazbPUOauujjCyPuBOshScxOr6brU+u4PKbv+p1iSIiJ00t/HbaroBVUx8lizsv/pIVW3m3cj4AFftfpOVIi7dFioicAgV+Ox1dASueTPPMzgy1/goGO1E2PvkHj6oTETl1Cvx2Or0CVkOU+5vnsyE+icibz5FMJvu4MhGRnlHgt9PZFbAAGmOwvOVidiXK2fTMH/uwKhGRnlPgt9PRFbDyJQnwRHQWma2rdEUsERlQFPjt5F8BqzP1mTJGZ2t5dc1zfViZiEjPKPA7UDW7kl9858pOQ39IyG3Zx/76KJmMWvkiMjAo8LvQWfdOU8LPHQ2fpjoa4Y0XX/CgMhGRk6fA70JX3TsNmTKWt1zMi8+9QDarqZNFpP/r1idtjTFXAP8GTAEOAf9urf2pMSYELAE+A6SBH1lrF+ftdx1wDzAWWAvcZK091LtP4fRqu8D5orueoabdkM0kAV5oOZOq9X9hygUXeVShiEj3nLCFb4ypBFYCdwHDgM8Bi40xfwvcCRhgMnA+cKMxZmFuv+nAg8BNwAjgTWB57z+FvtHZ+Pz6TBm1a/9LrXwR6fe606UzCXjYWvuotTZjrV0PrAEuAW4E7rbW1ltrdwM/AL6c2+8GYJW1dp21NgbcBlxijJnSy8+hT3Q2Pn+Yr5XRiWreWf9yH1ckInJyThj41toXrLVfabtvjBkOfATYjNtVsz1v8zeAGbnvp+evs9a2AtV56weUjk7ghoN+LhwdB+DwmofVyheRfu2kZss0xgwFHgdeATbmFrfmbdIKlOa+H9RuXfv1A0rV7EqAY2bRPH/qKF7ZFuCPTV+g3NfCVf/1LH//91d6XKmISMe6HfjGmA8Aj+G22q8H2vo48vs6SoEjue9b2q1rv37AaTuBC+/PqulOtOZQnxnEI+ubGX32e3x0zgRvCxUR6UC3hmUaYy7FbdX/AfiMtTZmra0HDuCetG0zlfe7cbbnrzPGlAITOLYLaMDqaFbNJAF++YfNHlUkItK1E7bwjTGTgSeAb1trH2i3+iHgDmPMq7hdON8E7sutexhYZ4ypAl4GFgObrbU7e6l2T3U6aieaZdH3nqG24f2Lp7S9KxAR8VJ3unS+BgzGHYq5OG/5T4DbgR8C23DfLfwMWApgrX3NGLMod38c7juEBb1XurdGlkeOG5ffpqbBXd528RRAoS8injth4FtrbwVu7WKTr+W+Otp3Je4Y/oKzcP60vD78NlnAOWa7eDLNstU7FPgi4jlNrXCK8qddcCA3/YLT4baddf+IiPQlXcS8B/JH7QAdTr8AXV9URUSkr6iF34s6+3DWwvnTPKpIROR9auH3orbW/s8f3UpTNEW5r4XPzZuq/nsR6RcU+L2sanYll543nt9//3vMSm+lecckspfPwnE67t8XEekr6tI5DXw+h8kfX0hzpoTBzbs5tPFPXpckIqLAP13OmzGJ14dVAXD4v5eRiWukjoh4S4F/Gs397Gd5NzWSkvQRdq3+jdfliEiRU+CfRuNHDaHG/B2ZLGS2PUO8Zo/XJYlIEVPgn2Yf/+RlbMoY/GR4e+VPNGe+iHhGgX+alUWClFd9npZMiJLDO6nbrBO4IuINBX4fuOzD5/BSyVwAap75FemWRo8rEpFipMDvAy9s3sOLLZP433Vf4J7DH2PlL//T65JEpAgp8E+ztitj1TfHabsy1vJdFfzxqRe8Lk1EiowC/zTr7MpYv1mzh0wi5lFVIlKMFPinWWdTIzekI+x/+ld9W4yIFDUF/mnW2dTI5b4W4q89S+s7W/u4IhEpVgr806yjKZNDQR8TgvUA7Hn0ATKxFi9KE5Eio9kyT7O2qZGXrd5Bbf37FzaPx6fz7rOvMjF2mD/8+j95ombCMes1pbKI9DYFfh9of2UsgGw2yw83f4L9B/7M73aWk0QXPheR00tdOh5xHIcvfG4ej0YvJNnudbftwuciIr1Jge+h0cNLiWWDHa7Thc9FpLcp8D1WMazjUTy68LmI9DYFvscWXjWNUPDYX0Mo4OjC5yLS6xT4HquaXcnXF3yIoYNCQJZy3xGuG7yBj0wd4nVpIlJgFPj9QNXsSn5z53yu+fAkFpa9wPnONvat/CHZTPrEO4uIdNNJDcs0xlwAPGGtHZW7HwaagUTeZi9Za6/Mrb8OuAcYC6wFbrLWHuqNwgvRjVfP4I63P87I6HKGVG+j7vnfMmLeQq/LEpEC0a3AN8Y4wM3AD9qtmgHUWWvHdLDPdOBBYD6wAbgXWA5c1pOCC1ko6OfrN36U++8/yM3hp2n8y2OEx01h0NSLvC5NRApAd7t07gS+CtzVbvlsYEsn+9wArLLWrrPWxoDbgEuMMVNOqdIicUbFIK5ZcBWPt84G4OBjD5Coec/jqkSkEHQ38Jdaa2fjttTzzQJGGWNeNcYcNMasMMaMy62bDmxv29Ba2wpU474rkC5cfO4ZDLvwajbGJ+Gk4uxbfo+ukiUiPdatwLfW7utkVQvwIjAPMEAUeDS3bhDQ2m77VqD05MssPjdefQ6bK67m3dQIMk017F9xL9lU0uuyRGQA69EoHWvtrdbar1tra6y1DcCtwPnGmErcF4P2nx4qBY705GcWi4Dfx7duvIjf++bTkCklsddS89T/I5vNel2aiAxQPQp8Y8x3jTH5nxAK5W5juN05Jm/bUmACed080rXyISX8w//4KL+KXU48G+DIa2tpeHGl12WJyADV09kyzwXmGGM+n7t/H/CktbbGGPMwsM4YUwW8DCwGNltrd/bwZxaVyeOH8dnrruD+36RooYT6RwOMfOZxbvzkeZpNU0ROSk8/eHUzUA+8BezGHY//BQBr7WvAImApUAucAyzo4c8rSslUmv1UUJ8ZBDjUtmR54JHNrNlY7XVpIjKAnFQL31q7BhiWd/8wcH0X268E1AfRQ8tW7yCdObbvPpHK8utVr6qVLyLdpqkVBoDOpkqubU4SP7i7b4sRkQFLgT8AdHUh9P2/vZNE7Z4+rkhEBiIF/gDQ0YXQAWYFd5GJNrHvN/9Ksm5/3xcmIgOKAn8AqJpdyS0LZlJRHsEBKsojLLpmOq+FP8SbydFkWurZ99t/JVl/wOtSRaQf00XMB4iOLoQ+24zm9p9kuIGnOKuphn0P/Qtjr/9XQiPGdfIoIlLM1MIfwCaMGcJ3vnwpy5If463kKNLNdex76F9IHNJkayJyPAX+AHf2+GF8+0uX8lDyb7HJsWRaGtn3m9uJ73/b69JEpJ9R4BcAM3E4t39lLg+nrmRbYhyZaDP7Hrqd1ne2el2aiPQjCvwCMaWynDu/OpdHMleyIX4m2WSMA/91N82vrfW6NBHpJxT4BeSscUO5639eypP+eTwXnQ6ZNDWP30/Dy3/QLJsiolE6hWbi2CHce8ul3PGzIE3NpVxbtoG6Pz3E2m11PLZ/HLUNMUaWR1g4f5qmZRApMmrhF6AxI8r4/tc/wt5RF/OL5rn8JTaZX+8YSk1DjCxQUx9lyYqtmnxNpMgo8AvU0EFh7v7qJQQmn88jrX9Dst2buXgyzbLVOzyqTkS8oMAvYJFwgH9ZdCFpjp+WATqflE1ECpMCv8D5/T4qOpl8bXg4RSYR6+OKRMQrCvwi0NHka0FSXBV4ib2/+JY+pCVSJBT4RSB/8jVwL5AeIc64QAPJw/vY+6vbaHjp92QzaY8rFZHTScMyi0T+5GvJVIZfPrGNH74Q5hOlm5hb8gZ1z/+W1rc2UXH11wgOH+txtSJyOqiFX4SCAR9f+tQMvnXTRTzLJSxtnkdzNkKsegd7fn4rDX95XK19kQKkwC9iF804gwe++VFCk2ZyT8MnWB8/i2wqQd1zv2bfr79N/MA7XpcoIr1IgV/kRgyN8N0vXcznPjGHRxJz+WnzZTRmy4jve5O9v/gnap/+OenoEa/LFJFeoMAXfD6HT82dzAPfrMI3YSZ313+C52PTyGSzNG18muqlX6dp87Pq5hEZ4BT4ctQZIwdx91cuYdGn5/Bs5iK+33g1b6dGk2ltovappez5+a20vrlRE7GJDFAapSPH8Pkcrrr4TC6aMZZfPbGd+zcM47zQbj5ZtoXy2j0ceOQeSiacQ/nczxKZcI7X5YrISVDgS4fKB5fwjc/N4soLJ/KzPwzje3sn8AH/fvZny2moK6X81S1cW/nfzL/2Skoqp3ldroh0g7p0pEvnnDWC//sPc7nyorPYkR5PQ6YMcKjPDOKhd8fz+M//g33LvkPLzvVksxmvyxWRLqiFLyfk8zlseOPQccuTBFgVnc2c6pXEqncQHDmeoRdczaBzPoIvVOJBpSLSlZMKfGPMBcAT1tpRufshYAnwGSAN/Mhauzhv++uAe4CxwFrgJmvt8ckh/V5nM2s2ZEp5tGUO88p2MKR2D7VPLaXuTw8xeOZlDJn1t/rUrkg/0q3AN8Y4wM3AD9qtuhMwwGRgKPC0MWavtXaZMWY68CAwH9gA3AssBy7rpdqlD40sj1DTQeiXRUK8Hp7Nn+umcl5oN1WlO5kQO0TjK6tofGUVJRPOYfCHLqNs6kX4gmEPKheRNt3tw78T+CpwV7vlNwJ3W2vrrbW7cV8QvpxbdwOwylq7zlobA24DLjHGTOl52dLXOppxMxz085VrZ/Dz/3M53/j8HOorzuOHDR/jB41X8Up8MikCxN7bRs3jD/DufV/k0Kqf0Lprq8bzi3iku106S621txtjqtoWGGOG4XbVbM/b7g1gRu776bgtewCsta3GmOrc+jd7UrT0vbaJ15at3kFtffS46+JWza5k7qzx2HfrefKlXTyyZRS/bz2f80K7uSTyNpXxGo68+ieOvPon/GXDKJt2MWVTL6SkchqOr+MLtIhI7+pW4Ftr93WweFDutjVvWStQmre+lWPlr5cBJn/GzY44jsPUScOZOmk4X7zmgzy/cQ/PrR/JD/Z/gNG+BmaFd3NBZDfDWxpo2vAUTRuewlc6hLIp51N69mwiZ56LL9zxxVpEpOd6MkqnJXeb/z+0FDiSt779/9789VLAhg4K86m5k/nkpWfx9t5Gnt9Qzbqto1ldN5MJ/sPMDL3LeSXVjGhtonnrczRvfQ78ASITziFy1kwik84lNHoijqORwyK95ZQD31pbb4w5gHvSdm9u8VTe7+LZnlsHgDGmFJjAsV1AUuAcx+Hs8cM4e/wwFl3zQbbvOswLm/fy4rZxrKqPMdbfwAeDe5gR3sMEaoju2kp011YAfKVDiEw8h5LK6UQmnkOwolIvACI90NNx+A8BdxhjXsXtwvkmcF9u3cPAuly//8vAYmCztXZnD39ml9ZsrO60n1m85fc5HG6IsuGNg9Q3xRk2uITKyqm8WT+RZ/c3UebEmBrchwnuZ1roAENam2jZ8TItO14GwFdSRviMKYTHfYCSM6YQHjsZf9lQj5+VyMDR08C/HfghsA13xM/PgKUA1trXjDGLcvfHAa8AC3r487q0ZmM1S1ZsJZ50R4HU1EdZssJtLSr0vdf+99PQHGfrm7XcsmAmH5w8kk32EJvsIZ54q5aH6+OM8jUxOXiQswMHmRI6xNBYC9F3thB9Z8vRx/QPHkF4zFmERk8iNGoioVETCJaP0Ylg6bey2SxkUmQSceKxOIlYlHg0RjIW5aU36li1uYmG1gwVwyIsvKp3G6xOf5350BgzCdj13HPPMX78+G7ts+iuZzocK15RHuEX37mydwuUk9bd308mk+W9g828+mYN23fVsWP3Yeqa4gx1WpgUqGVSoIaJgVrGB+oJO8njHs/xBwmOGEtwxDj3a/hYguVjCZaPwVc6BMdxTuvzlIEnm0mTTSXIJhOkk3Hi0VYS0RiJWJRENEYyHiMVi5FKxEklYqQTcTKJOJlkgkwqTjaZgHQC0kmcdAInncTJJPFlkvgzKfzZFH5SBLIpAqTxOcfn7ob4JJa3XEwyrx0eDvq5ZcHMkwr9PXv2MG/ePIAzc8PljyqoqRU6+zRoZ8ulb3X39+PzOUwaO4RJY4dwzaWTyWazHKxrxb5bz87qena+18DqvY0kkykqfE2MD9Rxhr/e/Qo0UE4LiUPvkTj03nE/ywmVEBhaQXDoKAJDK/APHkFgyHACg0fgH1SOf1A5vnCpXhT6ifwgziTjZJNxMslcq7g1SjIWJRGLkozFSMVjpBMxUokEmUTM3TaVgGQC0nFIJY6GsBvESfzZFAGSBOh6HigfEMp99UjuzyqV9ZHM+km6LwOknQCPReccE/YA8WSaZat39Forv6ACv7NPg44s11C//uBUfz+O4zBmRBljRpQxd5b7bi+dybKv5gi79jXyzt5G3jvYzPoDzRysayVMklH+Rkb7Gxntb2Kkr5mR/mYq/E1EEjGSNdUka6o7/3mBEP6yofhLh+ArHYK/dCj+yCB8kcH4SgbhLynDFy7FV1KKE4rgC0fwhSI4oRIcf7CgXyyy2SzZdJJsKukGcSrhfp90v29r7aaTCZK5IE7G46TisVzLOEE6EXNbxsk42VQCUm7r2Ekn8KWT+DIJfJkU/mySAN37kJ4/93WqMll3bqhk1k8iGyCFn5QTIO0ESTtBMr4gGV+ArD9E1h8EfwgnEIRAGF8whC8QwhcK4w+F8QfDBMMl+MPubTAcJlRSQigSIRSJEC6JUBIOEgr68fne/1tp/sfHOqytNxusBRX4C+dPO6aPGNy3RAvna/re/qA3fj8dnZS/6er35+WPxVPsq21hb80R9tUcYV9tCxvqWjlwuIW6higlJBjuO8JwXwvDfS0M9bUyzNfCMF8rQ3xRhviihFMJUo01pBprTv5JOg4EwoI4dg0AAAgUSURBVDjBME7ADQVfwA0EJxBwXxD8AfAHcHx+93uf3z3n4PO5t44Px+cDxweOk3sBcdzHxjnaSjwqm/snmyWbzZLJpMlmMmQyGbLpDJl0mkwm4y5Pp8ik3dtsOk0mnXK/z60jd59MCifz/q0vk8LJpvBnT/1T0qcSyvlBHM8Gjn6fJJAL4gBpX4iML0jWH4KA++W0BXGwBF8oTCAUwh+OEAiVEIyUEAyFCUdKj4ZwaaSEjW8c4nd/epPDjd4M+OiLBmtBBf6JPg0q3urp76c7J+VLwgHOGjeUs8YdP3onmUpT2xCjtiFKTUOUw41R6pvjvN0Uo64xRlNLnIYjCVLRVgb5YgxyYkdvS30Jypw4pU6CEidBxEkQ8SUoIUXYSea+Um7XQDJGNhmjrZe2P04kkXvpOKX50ZO57og0fhK58HW7KAIk8R/tqsi0tYz9QbK5QHbaAjkYwhd0W8P+UBh/yG0JB0rcAHZbxKWUlISJlAQZGvJTEgpQEvITDvrx+3t3eO6ajdX8ctU2Twd89EWDtaACH078aVDxVk9+P8tW7zjmPwOcXB9nMOBn7Mgyxo4s63K7ZCpNc2uS5pYETa0JjrQmaImmaI0laYkmaYqn2LyngZ3vNRBPpgkGfFQMixAO+UklU2TTCUjG8WWS7sm7dBKy7otBwEkTIIOfDD7HvfWTPfq9jyw+sjhO7pZsLpzdl4+2+wD5p/2yeVvgOEffHeD43HcOjh/H54Dv/XcUjt/PvliEnY0RommHSACmj3KYNDKIPxjEHwzhDwYJhEIEQmGCoTDBUJBwOEg46CcU9FOWC+BQ0E84F8rhkJ9QwHdSXVvvv3NrONoQ+Jvp47q9f0/19G+rN/RFg7XgAl8KV1+dlA8G/Awf4mf4kI7n9F+zsZrVL+0+GhDJVIbDjbETjqbIZLKkM1nSmQzZLPx50x6WP2upbYwxYmgJC+ZN4eJzz3h/hyy5Hhw3OB3HPZ/RdutzwOc4OD4Hn+Pg87nLuhu0azZW86cVW4mn3ecRTcHrNX6qqmb0aaOpPwyn7i8DPk53g1UfW5QBo7O+zL4+Kd9Va7ArPp9DMOCjJBTgldf38/PHXqe2MQbA4cYYv1y1na07aygfXOJ+DXFvhw0OM2xwmKGDwgwpCzG4NMSgSJDSkiAl4QDhoJ9gwIff55xUq/pUn0dv6w919Je/rdNNgS8DRmdTNJ9sH+eajdUsuusZrvnHx1h01zOs2dj5iJ2O9EZrsD+EXH9p1faHOnrrb6u/U5eODBi90cfZG90HvTGaoj+EXH8Zxtwf6iiWAR8KfBlQetrH2Rsn53pjNEV/CLn+Moy5v9RRDAM+1KUjRaU3WtZVsyu5ZcFMKsojOLhTQ5zsx9/7QxdCbzyPQqqjGKiFL0Wlt1rWPW0N9pcuhP7Squ0vdRQ6Bb4Ulf7SfQAKOel7CnwpKv2lZS3iBQW+FB21rKVY6aStiEiRUOCLiBQJBb6ISJFQ4IuIFIn+fNLWD3DgwAGv6xARGTDyMvO4683058AfC3D99dd7XYeIyEA0Fng7f0F/Dvz1wEeA/fTPiwaJiPRHftywX99+hZPNZo/fXERECo5O2oqIFAkFvohIkVDgi4gUCQW+iEiRUOCLiBQJBb6ISJFQ4IuIFIn+/MGrU2KMmQksBc4F3gEWWWuP+wCCdMwYswj4KRDPW/w14D+BJcBncD8I9yNr7eK+r3BgMMZcADxhrR2Vux+ii+NnjLkOuAf3AzNrgZustYf6vPB+qINjGQaagUTeZi9Za6/Mrdex7ERBBX7uP9VjwI+BS4FPA88YYyZaa5s8LW7gmAX80Fr7z/kLjTGLAQNMBoYCTxtj9lprl3lQY79ljHGAm4EftFt1J50cP2PMdOBBYD6wAbgXWA5c1meF90NdHMsZQJ21dkwH++hYdqHQunSqgKC19sfW2qS1djmwDfist2UNKLOBLR0svxG421pbb63djfuf8Mt9WdgAcSfwVeCudsu7On43AKusteustTHgNuASY8yUPqq5v+rsWHb2Nwo6ll0qtMCfDuxot+wN3BaBnIAxxo/bFfYFY8w+Y8xbxph/NsaU47493p63uY5rx5Zaa2fjti4BMMYMo+vjNz1/nbW2FahGx/e4Y5kzCxhljHnVGHPQGLPCGDMut07HsguFFviDgNZ2y1qBUg9qGYgqcP9z/Ro4E7e/+avA13Pr84+tjmsHrLX7Olg8KHfb2fHT320HOjmWAC3Ai8A83G6yKPBobp2OZRcKqg8f9w8h0m5ZKXDEg1oGHGvtAWBu3qItxpgHcPtD4dhjq+PafS25286On/5uT4K19tb8+8aYW4EaY0wlOpZdKrQW/nbcV/x8Uzn2rbR0whhzjjHmznaLQ0AMOMCxx1bHtZustfV0ffyO+bs1xpQCE9Dx7ZAx5rvGmGl5i0K52xg6ll0qtBb+84BjjPkG7hC4T+P2ST/a5V7SpgH4R2PMHtyRDucB/wu4Bffk9x3GmFdx3zZ/E7jPq0IHoIfo/Pg9DKwzxlQBLwOLgc3W2p1eFDoAnAvMMcZ8Pnf/PuBJa22NMUbHsgsF1cK31iZwux8+DdQB3wY+Za2t8bSwAcJauxe4Bnf0SBOwEvietfZ3wO3A67jBvz63bqlHpQ5EnR4/a+1rwKLc/VrgHGCBN2UOCDcD9cBbwG7c8fhfAB3LE9EFUEREikRBtfBFRKRzCnwRkSKhwBcRKRIKfBGRIqHAFxEpEgp8EZEiocAXESkSCnwRkSKhwBcRKRL/H5UezHjbUYNmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(results.G, '-')\n",
    "plot(results2.G, '-')\n",
    "plot(data.glucose, 'bo')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differences in `G` are less than 1%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0              0\n",
       "2      -0.127982\n",
       "4      -0.191302\n",
       "6       0.154946\n",
       "8        0.26436\n",
       "10      0.207439\n",
       "12      0.132947\n",
       "14      0.035413\n",
       "16    -0.0815439\n",
       "18     -0.222312\n",
       "20     -0.371789\n",
       "22     -0.522334\n",
       "24     -0.664226\n",
       "26     -0.801867\n",
       "28     -0.934267\n",
       "30      -1.05701\n",
       "32       -1.1642\n",
       "34      -1.25663\n",
       "36        -1.346\n",
       "38      -1.43135\n",
       "40       -1.5119\n",
       "42      -1.58699\n",
       "44      -1.65612\n",
       "46      -1.71812\n",
       "48      -1.77288\n",
       "50      -1.82036\n",
       "52      -1.86066\n",
       "54      -1.89392\n",
       "56      -1.91501\n",
       "58      -1.92508\n",
       "         ...    \n",
       "124    -0.661724\n",
       "126    -0.631811\n",
       "128    -0.601917\n",
       "130    -0.572193\n",
       "132     -0.54277\n",
       "134    -0.513756\n",
       "136    -0.485244\n",
       "138    -0.457312\n",
       "140    -0.430024\n",
       "142    -0.403431\n",
       "144    -0.377577\n",
       "146    -0.352881\n",
       "148    -0.329311\n",
       "150    -0.306834\n",
       "152    -0.285416\n",
       "154    -0.265021\n",
       "156    -0.245616\n",
       "158    -0.227165\n",
       "160    -0.209635\n",
       "162    -0.192991\n",
       "164    -0.177201\n",
       "166     -0.16262\n",
       "168    -0.149158\n",
       "170    -0.136734\n",
       "172     -0.12527\n",
       "174    -0.114695\n",
       "176    -0.104943\n",
       "178   -0.0959529\n",
       "180   -0.0876686\n",
       "182   -0.0800372\n",
       "Name: G, Length: 92, dtype: object"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff = results.G - results2.G\n",
    "percent_diff = diff / results2.G * 100\n",
    "percent_diff.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's find the parameters that yield the best fit for the data.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use these values as an initial estimate and iteratively improve them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>G0</th>\n",
       "      <td>290.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k1</th>\n",
       "      <td>0.03000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k2</th>\n",
       "      <td>0.02000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k3</th>\n",
       "      <td>0.00001</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "G0    290.00000\n",
       "k1      0.03000\n",
       "k2      0.02000\n",
       "k3      0.00001\n",
       "dtype: float64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = Params(G0 = 290,\n",
    "                k1 = 0.03,\n",
    "                k2 = 0.02,\n",
    "                k3 = 1e-05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_system` takes the parameters and actual data and returns a `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params, data):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "    \n",
    "    params: sequence of G0, k1, k2, k3\n",
    "    data: DataFrame with `glucose` and `insulin`\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    # params might be a Params object or an array,\n",
    "    # so we have to unpack it like this\n",
    "    G0, k1, k2, k3 = params\n",
    "    \n",
    "    Gb = data.glucose[0]\n",
    "    Ib = data.insulin[0]\n",
    "    I = interpolate(data.insulin)\n",
    "    \n",
    "    t_0 = get_first_label(data)\n",
    "    t_end = get_last_label(data)\n",
    "\n",
    "    init = State(G=G0, X=0)\n",
    "    \n",
    "    return System(G0=G0, k1=k1, k2=k2, k3=k3,\n",
    "                  init=init, Gb=Gb, Ib=Ib, I=I,\n",
    "                  t_0=t_0, t_end=t_end, dt=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>G0</th>\n",
       "      <td>290</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k1</th>\n",
       "      <td>0.03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k2</th>\n",
       "      <td>0.02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k3</th>\n",
       "      <td>1e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>G    290.0\n",
       "X      0.0\n",
       "dtype: float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gb</th>\n",
       "      <td>92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ib</th>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>I</th>\n",
       "      <td>&lt;function interpolate.&lt;locals&gt;.wrapper at 0x7f...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "G0                                                     290\n",
       "k1                                                    0.03\n",
       "k2                                                    0.02\n",
       "k3                                                   1e-05\n",
       "init                  G    290.0\n",
       "X      0.0\n",
       "dtype: float64\n",
       "Gb                                                      92\n",
       "Ib                                                      11\n",
       "I        <function interpolate.<locals>.wrapper at 0x7f...\n",
       "t_0                                                      0\n",
       "t_end                                                  182\n",
       "dt                                                       2\n",
       "dtype: object"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`error_func` takes the parameters and actual data, makes a `System` object, and runs `odeint`, then compares the results to the data.  It returns an array of errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(params, data)\n",
    "results, details = run_ode_solver(system, slope_func, t_eval=data.index)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_func(params, data):\n",
    "    \"\"\"Computes an array of errors to be minimized.\n",
    "    \n",
    "    params: sequence of parameters\n",
    "    data: DataFrame of values to be matched\n",
    "    \n",
    "    returns: array of errors\n",
    "    \"\"\"\n",
    "    print(params)\n",
    "    \n",
    "    # make a System with the given parameters\n",
    "    system = make_system(params, data)\n",
    "    \n",
    "    # solve the ODE\n",
    "    results, details = run_ode_solver(system, slope_func, t_eval=data.index)\n",
    "    \n",
    "    # compute the difference between the model\n",
    "    # results and actual data\n",
    "    errors = results.G - data.glucose\n",
    "    return errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we call `error_func`, we provide a sequence of parameters as a single object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how that works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "G0    290.00000\n",
      "k1      0.03000\n",
      "k2      0.02000\n",
      "k3      0.00001\n",
      "dtype: float64\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0      198.000000\n",
       "2      -71.523600\n",
       "4      -19.535535\n",
       "6        4.898956\n",
       "8        4.423982\n",
       "10      17.420932\n",
       "12      11.905305\n",
       "14       7.909446\n",
       "16       7.454111\n",
       "18            NaN\n",
       "19            NaN\n",
       "20            NaN\n",
       "22       6.458757\n",
       "24            NaN\n",
       "26            NaN\n",
       "27            NaN\n",
       "28            NaN\n",
       "30            NaN\n",
       "32       4.670366\n",
       "34            NaN\n",
       "36            NaN\n",
       "38            NaN\n",
       "40            NaN\n",
       "42       0.314528\n",
       "44            NaN\n",
       "46            NaN\n",
       "48            NaN\n",
       "50            NaN\n",
       "52       4.166524\n",
       "54            NaN\n",
       "          ...    \n",
       "124           NaN\n",
       "126           NaN\n",
       "128           NaN\n",
       "130           NaN\n",
       "132           NaN\n",
       "134           NaN\n",
       "136           NaN\n",
       "138           NaN\n",
       "140           NaN\n",
       "142      6.365824\n",
       "144           NaN\n",
       "146           NaN\n",
       "148           NaN\n",
       "150           NaN\n",
       "152           NaN\n",
       "154           NaN\n",
       "156           NaN\n",
       "158           NaN\n",
       "160           NaN\n",
       "162      5.142209\n",
       "164           NaN\n",
       "166           NaN\n",
       "168           NaN\n",
       "170           NaN\n",
       "172           NaN\n",
       "174           NaN\n",
       "176           NaN\n",
       "178           NaN\n",
       "180           NaN\n",
       "182      1.798265\n",
       "Length: 94, dtype: float64"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error_func(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`leastsq` is a wrapper for `scipy.optimize.leastsq`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we call it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.9e+02 3.0e-02 2.0e-02 1.0e-05]\n",
      "[2.9e+02 3.0e-02 2.0e-02 1.0e-05]\n",
      "[2.9e+02 3.0e-02 2.0e-02 1.0e-05]\n",
      "[2.90000004e+02 3.00000000e-02 2.00000000e-02 1.00000000e-05]\n",
      "[2.90000000e+02 3.00000004e-02 2.00000000e-02 1.00000000e-05]\n",
      "[2.90000000e+02 3.00000000e-02 2.00000003e-02 1.00000000e-05]\n",
      "[2.90000000e+02 3.00000000e-02 2.00000000e-02 1.00000001e-05]\n"
     ]
    }
   ],
   "source": [
    "best_params, fit_details = leastsq(error_func, params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first return value is a `Params` object with the best parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>G0</th>\n",
       "      <td>290.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k1</th>\n",
       "      <td>0.03000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k2</th>\n",
       "      <td>0.02000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k3</th>\n",
       "      <td>0.00001</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "G0    290.00000\n",
       "k1      0.03000\n",
       "k2      0.02000\n",
       "k3      0.00001\n",
       "dtype: float64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second return value is a `ModSimSeries` object with information about the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fvec</th>\n",
       "      <td>[198.0, -71.52359999999999, -19.53553501808568...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nfev</th>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fjac</th>\n",
       "      <td>[[nan, nan, nan, nan, nan, nan, nan, nan, nan,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ipvt</th>\n",
       "      <td>[1, 2, 3, 4]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>qtf</th>\n",
       "      <td>[nan, nan, nan, nan]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>cov_x</th>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mesg</th>\n",
       "      <td>The cosine of the angle between func(x) and an...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ier</th>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "fvec     [198.0, -71.52359999999999, -19.53553501808568...\n",
       "nfev                                                     5\n",
       "fjac     [[nan, nan, nan, nan, nan, nan, nan, nan, nan,...\n",
       "ipvt                                          [1, 2, 3, 4]\n",
       "qtf                                   [nan, nan, nan, nan]\n",
       "cov_x                                                 None\n",
       "mesg     The cosine of the angle between func(x) and an...\n",
       "ier                                                      4\n",
       "dtype: object"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit_details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have `best_params`, we can use it to make a `System` object and run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'The solver successfully reached the end of the integration interval.'"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(best_params, data)\n",
    "results, details = run_ode_solver(system, slope_func, t_eval=data.index)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results, along with the data.  The first few points of the model don't fit the data, but we don't expect them to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap08-fig04.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhU5fXA8W8yWUgIhkDCvm9HQNlBRNEU3ELFahF3EdFWK7a21tZ9wapo1bqAViwqYkV+4i6IqFQUFBQiixo87BCWkEDCmpD998edhMkyYQLJzCQ5n+fJMzN3mXtyCXPmfe973xNSXFyMMcYYE2xCAx2AMcYYUxlLUMYYY4KSJShjjDFByRKUMcaYoBQW6ABqm4hEAoOBXUBhgMMxxhhTlgtoDSxX1VzPFfU+QeEkp8WBDsIYY0yVhgNLPBc0hAS1C+DNN9+kVatWgY7FGGOMh7S0NK6++mpwf1Z7aggJqhCgVatWtGvXLtCxGGOMqVyFSzA2SMIYY0xQsgRljDEmKFmCMsYYE5QsQRljjAlKfh0kISIXAo8BnYF04J+qOs19r9JBIM9j829V9Tz3fpe592sNfAWMV9V0f8ZujDHGv/yWoESkNfAOcImqzheRAcA3IrIcpyWXqaoVxoGLSC/gFSAJWAE8AcwGRvgr9hKLklOZOX8te7JyiI+LYlxSTxIHtvd3GMYY0yD4rYtPVXcBCe7kFAo0BwpwWk4DgVVedr0G+FhVl6jqEeBu4AwR6e6PuEssSk5l6pzVZGTlUAxkZOUwdc5qFiWn+jMMY4xpMPx6DUpVD4pINJALfAa8oKrrgQFACxFZIyK7RWSOiLR179YLSPF4j2wgFTjVn7HPnL+W3Pyyw/Rz8wuZOX+tP8MwxpgGIxCDJI4AjXGmIJogIjcAh4FvgJGAADnA++7tY4Dscu+RDUT7JVq3PVk51VpujDHmxPh9JglVLcIZDLFCRF4GfqOqF3luIyK3Axki0h4neUWVe5to4JA/4i0RHxdFRiXJKD6ufGjGGONfqamptG9f/66H+60FJSJni0hyucWRwD4ReVhEenosj3A/HsHp3hOP94kGOuDR7ecP45J6EhnuKrMsMtzFuKSeXvYwxjQEL730ErfffnuNvuf27dsREQ4cOHDMbd98800ef/zx0tf9+/dHVWs0nkDxZwtqFdDW3Tp6DjgNuAG4BPgLMEhErnJv+xwwT1UzRGQWsEREEoGlwGRgpaqu82PspaP1bBSfMcbTzTffHNDjZ2ZmUlxcXPp65cqVAYymZvlzFN9+YBTwWyATeBm4UVW/wklUWcAGYAtOF+C17v1+BCYALwF7gN7AWH/F7SlxYHteve88Pnr6N7x633mWnIxpQIqKinjssccYNmwYp59+OjfccAPbtm1jypQp3HLLLQBMmTKFu+++m1tuuYX+/fszevRoVq1axZ/+9KfS1yWtG8/9oOpW04IFCxgzZgxDhgxh8ODB3H333eTn57NgwQKmTZvGokWLuOgi50qJiLB2rTN466effuLaa69l0KBBnH/++bz55pul73nttdfyzDPPcMkllzBgwACuuuoqNm7cWGvn73j49RqUqv4AnFnJ8r3A1VXs9y7wbi2GZowJEpOmL2PF2t1+Odagni158MahPm37+eef8/XXXzN//nyio6N54IEHmDZtWoUyPh999BEvvPACU6ZMYeLEiVxzzTW88MILPPXUU9x1111MnTqVKVOm+Bzjjh07+Pvf/84rr7zCoEGD2Lp1K5dffjlffPEFSUlJrFu3jrVr1/Liiy+W2S8zM5Px48czceJEXn31VdatW8dNN91EbGwsF154IQAffPABM2bMICEhgdtuu40pU6bw7LPP+hxbbbOpjowxxgdNmjQhPT2dDz/8kLS0NB599FEeffTRCtv16dOHxMREXC4XQ4YMoWvXrpx99tlEREQwbNgwtm/fXq3jJiQkMHfuXAYNGsTBgwfJzMwkLi6O9PSqJ9NZuHAhCQkJXH/99YSHh9O7d2/GjRvHu+8e/a5/0UUX0blzZ2JiYjj//PPZtm1btWKrbQ2hHpQxpg7xtUXjb8OGDeOBBx5g9uzZPP7447Rv354777yzwnZNmzYtfR4aGspJJ51U5nVRUVG1jhseHs67777LO++8Q6NGjejVqxe5ubllrjtVJjMzkzZt2pRZ1rZtW3btOloXsHnz5qXPw8LCKCysUJIpoCxBGWOMD1JTU+nVqxezZ8/m4MGDzJo1iz//+c9cf/31ZbYLCQnx6f1CQ0PJz88vfb1v375Kt5s3bx4ff/wx7777Li1btgQovd5UldatW7Nz584Kv0N8fLxP8QUD6+IzxhgfLFu2jFtvvZWdO3cSExNDbGwsTZo0weVyHXvnSnTu3JmVK1eyadMmsrOzee211yrd7uDBg7hcLiIiIsjPz+eNN95AVUuTW0REBAcPHqyw39lnn01WVhYzZswgPz+flJQU3njjDUaPHn1c8QaCJShjjPHBmDFjOOecc7jssssYMGAAb7/9Ns8//7zPLabyzjnnHH79619z5ZVXMmrUKE4//fRKt7vkkkvo1asX55xzDmeddRbLli3jwgsvZP369QAkJiaydetWzj777DL7xcbGMn36dL744guGDh3Krbfeyo033sjll19+XPEGQsix+jHrOhHpBGxeuHAh7dq1C3Q4xhhjPGzfvp2RI0cCdFbVLZ7rrAVljDEmKFmCMsYYE5QsQRljjAlKlqCMMcYEJUtQxhhjgpIlKGOMMUHJEpQxxpigZAnKGGNMULIEZYwxxicHDhzwqcpvTbEEZYwxNWDEiBF88cUXgQ6j2u66665Ky4ZU5vzzz2fHjh21HNFRNpu5MaZOW5Scysz5a9mTlUN8XBTjknpatetakpmZ6dfjWQvKGFNnLUpOZeqc1WRk5VAMZGTlMHXOahYlp9bK8T777DPOP/98TjvtNO655x6uuOIK3nvvvQrblW9NeZZ3Ly4u5j//+Q+JiYkMHDiQG2+8kbS0NAC2bt3KzTffzJAhQxgxYgRTp06loKAAgNWrVzNmzBgGDRrEBRdcwPTp00vfPy0tjYkTJ3LaaadxzjnnMGPGDK+/Q0pKCmPHjqVfv35MmDChTNI5cOAAf//73xkxYgR9+/Zl9OjRLF26FIDf/va3AFxxxRV88skn5OXl8cgjj3DeeefRr18/zj33XObNm3ecZ7ZylqCMMXXWzPlryc0vW2QvN7+QmfPX1vixNm/ezB133ME999zDkiVL6NChAytXrqz2+7z99tu8+eabvPzyyyxbtow2bdpwzz33kJeXx4QJE+jYsSOLFy9mxowZfPLJJ7zyyisA3HfffYwZM4YVK1bw3HPP8eKLL5KamkphYSE333wzrVu35uuvv2b69Om89dZbfPDBBxWOnZeXxx/+8AcSExNZvnw5119/Pd98803p+ieffJKcnBzmzZtHcnIyZ555Jo888ghAaSKePXs2o0aN4tVXX+Wnn35izpw5/PDDD4wbN44HHnigNKHWBOviM8bUWXuycqq1/ETMmzePYcOGlZa1uOmmm5g1a1a13+fjjz/mmmuuoUePHgD87W9/Y/v27SQnJ7Nv3z7uuOMOwsPD6dChAxMnTuS5557jpptuIiYmhkWLFtGhQwcGDx7MihUrCA0NZfXq1Wzbto133nmHsLAwOnXqxPXXX8/s2bO5+OKLyxw7OTmZ7Oxsbr75ZlwuF8OHDy9TpuO2224jIiKCiIgIdu7cyUknneS1tPwVV1zBZZddxkknncTu3buJiori0KFD5OTk0KRJk2qfl8pYgjLG1FnxcVFkVJKM4uOiavxY6enppRVtwamc26pVq2q/z549e8qUYm/SpAk9e/Zk7ty5JCQkEB4eXrrOs0T7s88+y7PPPstdd93FgQMHSEpK4v7772fHjh3k5OQwdOjQ0v2KiorKlJ73PHZ8fHyZIovt2rUrLR+fnp7OY489xvr16+nUqRPx8fFeS8sfOnSIhx9+mNWrV9O2bVs6d+4McMxS9NVhCcoYU2eNS+rJ1Dmry3TzRYa7GJfUs8aP1apVK1avXl36uri4mN27d1e6bWhoKHl5eaWvPcu5t2zZsvSaEzhJ4/XXXycxMZH09HTy8vKIiIgAjpZoLygoYNOmTTz88MOEh4eTkpLCHXfcwcyZMxkyZAjNmzdnyZIlpe+ZmZnJkSNHKsTVokUL0tPTKSgoICzM+fjfvXs3LVq0AOAvf/kLY8aMYebMmYSGhvL555/z3XffVfo7Pvjgg3Ts2JEXX3yRsLAwUlJSmDt37jHPY3XYNShjTJ2VOLA9t47tS0JcFCFAQlwUt47tWyuj+EaPHs2yZctYvHgxBQUFvP7662USjadOnToxf/58cnJyWLduHQsWLCjzPv/973/ZvHkzeXl5TJ06lZSUFPr06UPLli15+umnyc3NZdu2bfz73/9m9OjRuFwu7rnnHmbMmEFhYSGtWrUiNDSU2NhY+vTpQ0xMDC+++CJ5eXlkZmZyyy238Pzzz1eIa+DAgTRv3pznn3+evLw8li1bxsKFC0vXHzp0iEaNGhEaGsrWrVt58cUXS0vLA4SHh5eWlz948CCRkZGEhoaSnp7O008/DVBm+xNlCcoYU6clDmzPq/edx0dP/4ZX7zuv1oaYt2/fnsmTJ/Pggw8ybNgwNm7cSJs2bcp0yZW48847ycjIYNiwYTzwwAOMGTOmdN2YMWO44ooruOGGGxg2bBhpaWlMnjyZ8PBwpk2bxubNmxk+fDhXXXUV5557Ln/84x8JCQnhueee44svvmDw4MGMGjWKoUOHcvnllxMREcHLL7/MmjVrGD58OKNGjaJbt2488MADFeIKCwtj2rRpJCcnM2TIEJ555pmSarYAPProo/z3v/+lf//+/P73v+eiiy4iPz+f1FRnVOSll17K7373O2bNmsW9997LkiVLGDhwIFdccQWDBw8mLi6OdevW1dg5t5Lvxhjjg507d5KdnU23bt1Klw0bNox//vOfnHnmmQGMrG6zku/GGHOC0tPTGTduHKmpqRQXFzN79mzy8vLo169foEOrt2yQhDHG+KBfv378/ve/59prr2X//v106dKFl156iZiYmECHVm/5lKBExAUMBAYBLYBCIA1Yrqqrai88Y4wJHuPHj2f8+PGBDqPBqDJBiUgc8CfgD0BzYBOwF3AB8UBHEdkFvAS8oKr7vL2X+/0uBB4DOgPpwD9VdZqIRABTgUtxkt+/VHWyx36XufdrDXwFjFfVyu8eM8YYUy94vQYlIuOAlUAH4AYgRlVFVYep6mmq2hVoBtwM9AZ+EpHrqni/1sA7wJ2q2gQYCzwrIgOASYAAXYHBwHXu4yMivYBXgPE4SXI9MPuEfmtjjDFBr6oW1GBgsKpmeNtAVQ8A84B57gR0L/C6l213iUiCqh4UkVCcZFMAHASuw2kVZQFZIvIUcBMwE7gG+FhVlwCIyN3ubbqr6vpq/r7GGGPqCK8JSlX/WJ03UtVdwK3H2OagiEQD+93HfgLIwOm6S/HY9BfgVPfzXsAKj/fIFpFU93pLUMYYU0+d0Cg+ERkMPKGqI6qx2xGgMdAH+AQomUgr22ObbCDa/Tym3Lry640xxtRDJzrMvBlw9jG38qCqRUAesEJEXsYZGQjgObtjNHDI/fxwuXXl1xtjjKmH/HajroicLSLJ5RZHAlk4Q9bFY/nJHO3yS/Fc5+4i7EDZLkFjjDH1jD9v1F0FtBWR24HngNNwRgdegpOgHhSRNThdene4twGYBSwRkURgKTAZWKmqNTfhkzHGmKDjtxaUqu4HRgG/BTKBl4EbVfUr4AHgJ+BnYDnwLs69Vajqj8AE9+s9OEPax/orbmOMMYHhtQUlIm/7sH+1qnWp6g9AhVkVVfUIMNH9U9l+7+IkLWOMMQ1EVV18h33Yf6P7p96as3Ad36zZyUM3nk7TJpGlyxclpzJz/lr2ZOUQHxfFuKSetTbNvzHGNERV3Qd1vT8DCVapuw+ycft+5i/dwpXnOWM1FiWnlqnimZGVw9Q5TqVNS1LGGFMzquriG+frm6jqzJoJJ/iMHNyBL5O38+nSLYwd2Z0wVygz568tU2IaIDe/kJnz11qCMsaYGlJVF9/T5V43A4qAnTj3MbXHmTR2Pc6URPVSn27xtE2IYUfGIb7/OY1hfdqwJyun0m29LTfGGFN9XkfxqWpCyQ/wELAI6KSqHVW1O9AO+BSY449AAyUkJIRRZ3QC4JNvNwMQH1f+vmGqXG6MMab6fB1m/iBwm6ruKFmgqnuAu3DKcdRrIwZ1IDLCxer1e0jdfZBxST2JDHeV2SYy3MW4pJ4BitAYY+qf6twH1bqSZV1x5tar12Kiwjm7fzsAPl26hcSB7bl1bF8S4qIIARLiorh1bF+7/mSMMTXI15kkXgNeF5GHcWpEheDMBHEv8GwtxRZUkoZ14rPvtrJw+TaudQ8pt4RkjDG1x9cEdTfOrOMPAi3dy3bhVMR9qjYCCzbd2jVFOsahW7P4auUOzh/aMdAhGWNMvVZVRd1h7sKCqGqRqj6kqq2BFkALVW3bUJJTiVHDOgPwyTebKS4uDnA0xhhTv1V1DeopIENEPhCRW0SkOziDI9wDJBqcM/u2ITYmgk0795OyOTPQ4RhjTL1W1TDzYUBnnBLupwDzRWSziLwsIpeKSJy/ggwWEeEuLhjaCYCPFtfrGZ6MMSbgqhzFp6oHVPV9Vb1FVbsBI4Bk4ApgvYh8JyKP+CPQYJE0rBOu0BCW/biL9MzyhX6NMcbUlGqV21DVzao6TVUvxbkW9ScawDBzT81jozizb1uKimHeN5sDHY4xxtRbPo3iq2JevmKcltSZwHeqml9jkQWxi87qwlcrt7Pgu61ceZ7QKNKfdR+NMaZh8PWTdTxwFk5rqaSSbXcgGtgKxAFZInKeqm6o6SCDTY8OcaVDzr9MTiXJPbrPGGNMzfG1i28VsABor6oDVHUAzlx8HwJvA/HAfOD5WokyCF00vAsAHy/ZZEPOjTGmFviaoK4H/qaqWSUL3CXc7wNuUtVCnBklzqj5EIPTsD5taB7biNTdh1i1LiPQ4RhjTL3ja4LKBTpVsrwzUFIYKRIoqIGY6oQwVyhJwzoBNljCGGNqg6/XoF4GXnPPxbcCJ7ENBO4HpotIC+AZ4H+1EmWQOu+0jsz+TFmekkZGVg4JVm7DGGNqjE8tKFV9ACcB3Q0sBb4B/g487l42ENgP3Fo7YQanuCaNGHZqG4qKYcGyLYEOxxhj6hWfx0er6uPA4yLSHMhX1QMeq+e7fxqcUWd05utVO/jsu61cfq4QHlatW8uMMcZ44XOCEpFEoDfOtSZEpHSdqv6rpgOrK3p1bkaHVk3YlnaQZT/tYni/toEOyRhj6gVfb9R9FvgjsI2KM0cUAw02QYWEhDBqWGdeem8Nn3y72RKUMcbUEF9bUOOACar6em0GU1f9amA7Xp/3Mz9t3Mu2tAN0aHVSoEMyxpg6z9cLJtnA97UZSF0W3SicxAFOdd35324JbDDGGFNP+JqgHgGeEhGb08eLknuiFq5IJSe3wdwOZowxtcbXLr61wGPABs/BESVU1VWTQdVFndvE0qtzM1I2Z7Loh+0knd4p0CEZY0ydVp0bdZcBr+F09x0XETkX596p7kA68KSqThORSOAgkOex+beqep57v8twEmRr4CtgvKqmH28ctWFRcirbdx8CYNp7a2gU7uJXg9oHOCpjjKm7fE1Q7YEkVd10vAcSkfbAu8B1OJPMDgQWiMgWYC+QqaqtKtmvF/AKkIQzi8UTwGyc4olBYVFyKlPnrCY335n1qbComClzVhESAokDLUkZY8zx8PUa1Oc45TZORCdglrtCb5GqLgcW4UwwOxBnxvTKXAN8rKpLVPUIzswVZ4hI9xOMp8bMnL+2NDmVyC8oYub8tQGKyBhj6j5fW1DfAVNFZAywAShTmFBV/36sN1DVxcDiktci0gwYDrwBXAC0EJE1QEvga+DPqroD6IXTcip5n2wRSQVOBdb7GH+t2pOVU+nyDC/LjTHGHJuvLahzgeVADNAPGOzxM6i6BxWRWOAjnMT3IXAYZ36/kYAAOcD77s1jqHjdKxunWGJQiPcySWzjRlZp1xhjjpdPn6Cq+quaOqCI9MBJSinA1apaBNxebpvbgQz3davDQPkMEA0cqqmYTtS4pJ5lrkGVCAmBwsIiXC6bn88YY6rL6yeniDwoIj7XjxCRJu5yHFVtcxZOq+kD4FL3NSVE5GER6emxaYT78QhOIhOP94gGOriXB4XEge25dWxfEuKiCAESmkbRNCaCQzkFfJ+SFujwjDGmTqqqBbUf+FlE3gHeU9Vl5TcQkRCcLr5rgN/ilOSolIh0BeYC96rqlHKr+wCDROQq9+vngHmqmiEis4Al7slqlwKTgZWqus6XX9BfEge2LzNi74OvNvLKRz8xd8lmTj+1TQAjM8aYuslrglLVZ93J6e/AZyJSgHPD7h5wGgo4s5uHADOAM1R1WxXHmgg0ASaLyGSP5S8ANwDP4wzACAPmAb93x/GjiEwAXgLa4rTAxlb7N/WzRhHOvctrNuxh3KRPmXBhbxtybowx1RBSXFx8zI1EpDGQiDMcvCVQBKQBycCXqppbizGeEBHpBGxeuHAh7dq188sxy98XBRAZ7uLWsX0tSRljjIft27czcuRIgM6qusVzna+DJA7jtGrm1Xh09VBl90Xl5hcyc/5aS1DGGOMjG15WC7zdF+VtuTHGmIosQdUCb/dFxTf1eVCkMcY0eJagasG4pJ5Ehlec4P203hWmGjTGGOOFTXVQC0quM82cv5Y9WTk0jgrnUE4+G7bvC3BkxhhTd/iUoNz3O12MM4ovHGdoeSlf5uJraDzvizqSW8D4f3zGL1uzWLctix4d4gIcnTHGBD9fu/ieBebglLwYwgnOxdfQNIoM44KhHQF4f9GGAEdjjDF1g69dfGOBm1V1em0GU5+NHt6FD77ayLdrdpK29zCtmjcOdEjGGBPUfG1BReJUsjXHqXlsFGcPaEdRMXy0+LjrPhpjTIPha4J6FbhDRGxQxQm4+OyuAHz+3VYOZucdY2tjjGnYfE04XYALgctEZCtQ5tNVVYfUdGD1Uec2sfTvkcDKdRl8unQLY0f2CHRIxhgTtHxNUGvcP+YEXZLYjZXrMvh48SYuPrsr4WEV75cyxhjj+1x8k2o7kIaiX48EOrU+iS27DvDVD9s5Z0jHQIdkjDFByeeZJEQkSUQWi0imiOwTkWUickVtBlcfhYSEcEliNwDeW7SBoqJjzyZvjDENkU8JSkTG4VTB/RG4Fae20w/AayJyde2FVz8N79eWhLgoUncf4rufdwU6HGOMCUq+XoO6G7hdVV/wWPamiKwB7gLerPHI6rHwsFB+m9iNae//yNsL1zP0lNaEhIQce0djjGlAfE1QHYFPK1n+BVWUeTeORcmppfPyxcdFMS6pJ+ee1pH/+3wdG1L3sWpdBv2lRaDDNMaYoOLrNaj1wMhKlp8DVFXmvcErqa6bkZVDMZCRlcPUOatZumYnF53VBYA5C9cHNkhjjAlCvragngBeFZHewDL3stOB3wG31EZg9UVV1XWn/PVXvPu/9fy4cQ9rN2fSs3OzAEVpjDHBx6cWlKrOAm4AhgGvAFOBAcBYVX2t9sKr+6qqrts4Kpxfn+m0ot5euM6fYRljTNDzeeoiVX0TGwxRbfFxUWRUkqRCQkO46K8f0iy2Ea7QEFas3c2mHfvp0jY2AFEaY0zw8ZqgROSfwCRVPex+7pXVg/JuXFJPps5ZXaGbr+T+p737j+AKdUbwzf5cuWe8zRpljDFQdQtqME5xQnBqQHm7o9TuNK1C+eq6IaEhFW7OLXS/XvrjLmtFGWOMm9cEpaq/8nie6G07ERsffSye1XUv+uuHVW771me/cO/1p/kjLGOMCWq+ziRRKCIJlSzvAFhxo2qIj4uqcv2yn9LYsH2fn6IxxpjgVdU1qCuBS9wvQ4DpIpJbbrOOQGYtxVYvebsm5WnK/63kub/+yut6Y4xpCKpqQX0OHAIOu1/nuJ+X/BwCvgMurs0A65vEge25dWxfEqpoSW3aeYD1qVl+jMoYY4JPVdeg9gATAERkC/CUqh72tr0vRORc4HGgO5AOPKmq00QkAufeqkuBQuBfqjrZY7/LgMeA1jil58eravqJxBJIJdekLvrrh15HmMxaoDx441C/xmWMMcHE53pQItJCRPoBJRX2QoBIYKBnMvFGRNoD7wLXAR8CA4EF7uSXCAjQFYgFPhWRHao6U0R64dwcnASswJnVYjYwwsffMWh5u0cKYMXa3Ta7hDGmQfN1kMQNQCqwGPgSWOR+nA9c5uOxOgGzVPV9VS1S1eXu9zkDJ2k9qqpZqroFeAq4yb3fNcDHqrpEVY/gzKx+hoh09/G4QWtcUk8iw8tW1I0Md3Fa71YAzJj3M8XFNorfGNMw+TpZ7N04XXAJwF6gBzAUUGC6L2+gqotV9eaS1yLSDBgOrMTpukvx2PwX4FT3816e61Q1GydZnkod53k9KgRIiIvi1rF9+cuVA2gSHUHK5kyWr90d6DCNMSYgfE1Q7YEXVHUvTkLprarfA7dxHJPFikgs8BHOIItk9+Jsj02ygWj385hy68qvr9MSB7bn1fvO46Onf8Or950HwB+f/pKD2XkAvDhnVemNvMYY05D4mqD2AyXDztYBfd3PFafrzmci0gNnRvTdOIMiDrpXeQ5ri8YZJQjOiMHyQ94819cbnqU5Suw9kMtL764OYFTGGBMYviaoz4F/iUgn4FvgChHpCFyNk2h8IiJn4bSaPgAuVdUjqpoFpOEMkihxMke79VI814lINNCBsl2C9UJlpTkAPvt+G3lV3DdljDH1ka+zmf8FmAlcCLyEUwdqM1AA/N6XNxCRrsBc4F5VnVJu9RvAg+4S8jHAHcBz7nWzgCUikggsBSYDK1W13tWn8Faao6iomHGTFpCdk19akbdk6iRjjKmvfE1QpwCXqGrJJ+gI9/Dvfaq608f3mAg0ASaLiOew9BeAB4CngZ9xWnUv4yRCVPVHEZngft0WpwU21sdj1ilVDTs/nJMPHK3IC1iSMsbUa74mqP3VZeUAAB9gSURBVLdx7jtaU7JAVavVxaaqtwO3V7HJRPdPZfu+i3MPVb3myzRIcLQiryUoY0x95us1qA0cHRhhakllw8698dYdaIwx9YWvLaj1wAwRuRvYiDMvXylV9fVmXXMMnqU5ACY88lml3X7HmhXdGGPqOl9bUAU4gyS+A/ZQdtLYE5qfz1Ststkmwl2hjEvqGaCIjDHGP3xtQb0GLFXVfM+FIhIJjKrxqEwpz4q8JS2pxlFhnNG3bSDDMsaYWudrC+pLoGkly7vgDAM3tahkton3/zmadi1i2Hcojw++2hDosIwxplZVVbDwD8Ak98sQIEVEys+5E4Mz9ZHxgzBXKDddcir3T1vK7M/XcVb/drRsVi9mfDLGmAqq6uL7D871pVDgVeAfOFMelSjGmW5oYa1FZyro16MFZ/Vry9erdvDy+z9y/w2nBTokY4ypFVUVLCwZGIGIbAa+cS8zAXbDb05hxS+7+T4ljWU/7WLoKa0DHZIxxtQ4XwdJfA1cLCIDgXCcLr9Sqvr3mg7MeNfspEZcc0FPXv7gR6a9/yN9uycQFenrP6UxxtQNvg6SeBaYg1PVdggw2ONnUO2EZqoyalgnurSNZc++HGZ/poEOxxhjapyvX7vHAjerqk/FCU3tc7lCmXhpX+54/ms++Hojw/u1pVv7ygZaGmNM3eRrCyoS+Ko2AzHV16NDHKPP7EJRUTHP/d9K8guKAh2SMcbUGF8T1KvAHSJiFzqCyKLkVL79cRcAW3Yd4Kk3VwQ4ImOMqTm+JpwuOLWgLhORrUCe50pVHVLTgZmqlVTf9Zz5/Ns1u3hn4TouHdkjgJEZY0zN8DVBrcGj1IYJPG/Vd99c8AuXJHbD5fK1cWyMMcHJpwSlqpOOvZXxJ2/lNgoKi3lv0QbGWivKGFPH+XxNSUROA/4K9ABGA1cCm1V1Ti3FZqpQVfXdWQt+ob+0oFs7G9VnjKm7fOoHEpFRwP+ALEBwbtYNA/7rLsdu/KyyMhyR4S769YinoLCYp99M5kieTfxhjKm7fL1Q8TDwZ1W9Cac2FKr6GHArYLNIBEBl1XdvHduX+yYMpX3LGLanH+L1uSmBDtMYY46br118PYEvKlm+EHi+5sIx1VG++m6Jv141kDue/5q532wmIjyUxat3sicrh/i4KMYl9ax0H2OMCTa+tqC2U/mURucCW2suHFMTurZrytUXOBV331u0kYysHIqBjKwcps5ZzaLk1MAGaIwxPvA1QT0KTBORuwEXcKGI/At4BniytoIzx++SxG6Eh1X8583NL2Tm/LUBiMgYY6rHpwSlqjNxRu2dj1MjahIwFLhKVV+pvfDM8XKFhnid+sjbEHVjjAkm1Zm66DNgharuBRCRwVg13aCW4GUoenxcVACiMcaY6vF1mHkvYCNwl8fij4EfRaRrbQRmTlxlQ9EjwkMZl9QzQBEZY4zvfL0GNRVYjNO1V6Iz8J17nQlCpUPRmx5tMcXHRnFG37YBjMoYY3zja4IaDDysqodKFqhqDvAIcEZtBGZqRuLA9rx6/3nMfOh8msc2Yueew/znwx8DHZYxxhyTrwkqEzilkuU9gIM1F46pLXFNGnHP+CGEh4Uy/9stfLp0S6BDMsaYKvk6SGI68LKItAdWAMXAAOB+oNqj+ERkCDBXVVu4X0fiJDrPMh7fqup57vWXAY8BrXEKJ45X1fTqHreh69EhjomX9uXZ2SuZ9v4aOrRqQq/OzQMdljHGVMrXBPWoe9v7gQT3snSqeR+UiIQANwBPlVt1KpCpqq0q2acXThJMwkmOTwCzgRG+HtccNXJwBzbt2M9HizcxecZynv7zWbSIiw50WMYYU4Gv90EVqeqDqtoSaAE0VdVWqvqEqlanzvgk4A841648DQRWednnGuBjVV2iqkeAu4EzRKR7NY5rPEwY3Zu+3ePZdyiXh6cv43BOfqBDMsaYCnyuaici/UXkRmAscI2I3FLyU43jvaSqA3FaQp4GAC1EZI2I7BaROSJSMtSsF1A666mqZgOpOK0ucxxcrlDuGjeYdi1i2Jp2kMdfX05BYXW+ZxhjTO3z9T6oe4Fk4HGc2cv/5vFzh68HU9WdXlYdBr4BRuKU88gB3neviwGyy22fDVi/1AmIiY7gwRuH0jQmklXrM3jxndUUFxcHOixjjCnl6zWoG4H7VfXR2ghCVW/3fC0itwMZ7kEZh4HyUx9EA4cwJ6RV88bcf8Np3P3iN3z+/TZaNovm8nMl0GEZYwzgexdfc+Dt2gpCRB4WEc/pDSLcj0dwuvfEY9tooAMe3X7m+PXoEMcdVw8kJAT+++kvzF+6JdAhGWMM4HuCege4uhbj6AM8LSJNRaQp8BwwT1UzgFnAb0Qk0T0cfTKwUlXX1WI8DUpuXgGNG4UD8OI7q5n23poAR2SMMb538eUAd4vIWGA9Ze9XQlUvO8E4bsApfLjBHdM84Pfu9/7RXVb+JaAtzvRKY0/weMZtUXIqU+esJje/sHTZ3G82ExHh4voLewcwMmNMQ+drgorGacnUCFVdBDT1eL2XKlpoqvou8G5NHd8cNXP+2jLJqcR7X25gSK9W9O5iN/IaYwLDpwSlqtfXdiAmMKqqDTVp+lIe+t3pNtuEMSYgfK4HJSKn4gwx741z7eoX4HlV/baWYjN+EO+lZlRkuIuc3EIe+s9SJv1uGD07NwtAdMaYhszX+6CSgB9wRvO9A8wBmgBfich5tReeqW2V1YyKDHdxy5g+nNW/LTm5hTz4n6X8siUzQBEaYxqq6szF94iqetaDQkTuA/6BU23X1EGJA9sDzrWoPVk5xMdFMS6pJ4kD23P2gHZQDF+v2sEDLy/l/htO49Su8QGO2BjTUPiaoHoClY3Umw3cU3PhmEBIHNi+NFF5crlCuf2qARACX6/cwUMvL+XO6wYzpFeFOX2NMabG+Xof1DagfyXLB+LMam7qKSdJDSTp9E7kFRTx6Gvfsyg5NdBhGWMaAF9bUC8AL4lIO2CZe9npwL3AP2sjMBM8XKEh/GFMH2Kiw5mzcD3/eusHDmbnM3p4l0CHZoypx3wdZv68iDTBKXVRchFiJ/Cgqk6treBM8AgJCWHcqF40bhTOjHkpvPzBj+zOzOb60b1xhYYEOjxjTD3k8zBz90Sxj4pICyBHVa3UewM0ZkR34k6KZMrbq/jw642kZ2Vz+1UDaBTh85+SMcb4pMprUCISKSI3i0hcyTJ3qfXfi8hEEYmoYndTT40Y1IGHfnc6jRuFsfTHXdz772/IPHAk0GEZY+oZr197RSQWZ/h4H5xqt8s8VrcBbgKuFJEka001PH27J/DPPw5n0vRlrNu2j7888xXnn96RL77fVmG4ujHGHI+qWlD348zB111VPZMTqvpX4BQgARtm3mB1aHUST912Fr06NyPzwBHeWqBkZOVQDGRk5TB1zmob8WeMOW5VJajfArer6vbKVqrqFuBO4NJaiMvUEXFNGvHIzWfQKNJVYV1ufiEz568NQFTGmPqgqgTVCjhWzaVVON19pgELDwslN7fijOhQ9WS0xhhTlaoSVCrQ4xj7dwfSai4cU1fFx0VVurxJYxtHY4w5PlUlqLeBh7yN1HMvfwiYWwtxmTqmsklnAQ4czuOZt37gcE5+AKIyxtRlVd28Mhm4GEgWkeeBFcB+IA4YDPwRcOFMFmsauAqTzjaNok/3eBav3MH/VqTy48Y9/PmK/vTplhDgSI0xdYXXBKWq2SIyDHgCeBKnvAZACLAXeAP4h6pm1XqUpk6obNLZMb/qzjNv/cD61H3c++9vGT28C9cm9SQq0m7sNcZUrcobdVX1oKreArTAKVR4JiBAC1W93ZKTOZb2LZvw5B+Hc/UFJ+MKDeHjxZuY+OT/WLF2d6BDM8YEOV/n4svDqaBrTLW5XKFcca4wuGdLpsxZxcbt+5k0fRln9m3D7y4+lWYnNQp0iMaYIORruQ1jTljXdk15+k9nccNFvYmMcLFk9U5ufnwh7325nvyCyoepG2MaLktQxq9crlAuPrsbL/xtBEN6tSInt4DX5qYw8ckv+T4ljeLi4kCHaIwJEpagTEC0bBbN/TecxqTfnU67FjHs2nOYf7zyHff8+xt+2ZoZ6PCMMUHAEpQJqAEnt2DKHb/ixt+cQpPocH7auJe/Pb+YR1/7jm1pBwIdnjEmgGysrwm4MFcovzmrKyMHd+C9L9fz4debWPZTGt/9nMawU9tw+bk96NwmNtBhGmP8zBKUCRoxUeGMG9WLZrGNeH1eCkdyC/lmzU6+WbOT03q3YsyvutOzc7NAh2mM8RNLUCaoLEpOZcbHKeTmlx3V993PTotKOsZx8dldOf2U1rhc1kNtTH1mCcoElZnz11ZITgDRjcIIDQlBt2bxxMwVJMRFccHQTpw7pANxdh+VMfVSQBKUiAwB5qpqC/frCGAqTm2pQuBfqjrZY/vLgMeA1sBXwHh36XlTz3grz5FzpIC3H/s1/0tO5cOvNrJzz2HemL+WWQt+4fRTW3P+0I6c2i0BV2iInyM2xtQWvyYoEQkBbgCeKrdqEs4USl2BWOBTEdmhqjNFpBfwCpCEM2HtE8BsYITfAjd+Ex8XRUYlSSo+LopGkWGMGtaZC4Z2YtW6DOYv3cz3P6exZPVOlqzeSXzTKH41sB0jB3egbUJMAKI3xtQkf7egJgG/Bh4B7vNYfh1OqygLyBKRp4CbgJnANcDHqroEQETudm/TXVXX+zV6U+vGJfVk6pzVZbr5IsNdjEvqWfo6NDSEASe3YMDJLdizL4fPv9vKwhWp7M7MZs7C9cxZuJ5u7WIZ3q8tZ/ZtS4tm0YH4VYwxJ8jfCeolVX1ARBJLFohIU5yuuxSP7X4BTnU/74XTcgJKZ1lPda+3BFXPVCjbERfFuKSeFWZJLxHfNIorzz+Zy88VUjbvZeHyVL5Zs5MN2/ezYft+XpubgnSI47RTWjH0lNa0b9mk0vcxxgQfvyYoVd1ZyeKSvphsj2XZQLTH+mzK8lxv6pnKynYcS2hoCKd0jeeUrvH8YUwfkn/Zzdcrd/B9ym50Wxa6LYuZn6ylbUIMg3u1ZODJLejdpTnhYRWLLBpjgkMwjOI77H70rBkeDRzyWF++nrjnemPKiAh3cfqpbTj91DYcyS1g5bp0lv2Uxvc/p7Ej4xA7vjrEB19tpFGEi1O7xdO3ewJ9usXTsdVJhNogC2OCRsATlKpmiUgaziCJHe7FJ3O0yy/FvQ4AEYkGOlC2S9CYSjWKDCtNVgWFRazdnEnyL7tJ/iWdLbsOsDxlN8tTnNpUsTER9O7SnF6dm9OrczO6tIm1e62MCaCAJyi3N4AHRWQNTpfeHcBz7nWzgCXu61ZLcUrRr1TVdbUd1KLkVJ+vhZjg4+3f79Ru8Yy/sDd79+ewal0GazbsYfX6DPbuP8K3a3bx7ZpdAERGuOjWrik9OsTRo0NTurVrSstm0YSEWCvLGH8IlgT1APA08DPOBLYvAy8BqOqPIjLB/bot8B0wtrYDWpScWmY0WUZWDlPnrAawJFUH+PLv1zw2ipGDOzBycAeKi4vZuecwKZv2krI5k5TNe9m55zA/b9rLz5v2lr5v46hwurSJpUvbWDq1bkKHVifRoWUTGlkJe1PPFRUVk5dfSF5BEXn5hSxevZ33F20k60AuCbX0BT6kvtffEZFOwOaFCxfSrl07n/eb8Mhnld6PkxAXxav3nVdzAZpaURP/fvsP5bI+dR/r3IMsNm3fz75DuZVu26JZNO0SYmjbIoa2CTG0iW9M6/jGJDSNsm5CU+OKi4vJLygiN7+QvPxC92ORx/OS5UWlzz2Ti7ft8/KLyCs4um++x/sUFBZVGVNkuItbx/atdpLavn07I0eOBOisqls819nXPi+8zWjgbbkJLjXx7xcbE8mgni0Z1LMl4HwoZB44wqYd+9m0cz/bdh1k2+6DbE8/SHpmNumZ2fxQboITV2gICXFRtIiLLn2MbxpF89hGxMc6j42jwq3bsB4oSRolH/q5ee5Hj+eliSPP+dDPzS/wWFfksa6wzPuUSTjuBBKItkVEuIuIsFCyj+RTVO74ufmFzJy/tkZbUZagvKhqRgMT/Grj3y8kJITmsVE0j41icK9WpcsLCovYteewM0Iw/RA7Mg6xa+9h0vYcZs/+I6TtzSZtb/k7JY4Kc4USd1IkcU0iiY2JJLZxJLExEZzUOIIm0RHEREfQJDqcmOgIohuF0bhROFGRYTbi0AeFJd1SFVoHR5d5JoHKEkNJ8vB8XT7plDz6M2mEuUKJDA8lItxFZITLSR7hLiLdSaTM65Ltwo9uFxEeSkRY2fURlW3rfq/wsNDSL1IX/fXDSmOq6S/wlqC88GVGAxO8/PnvF+YKpX3LJpXeBJyXX0h6VjbpWTlkuB/37Mth7/4c9u4/wt79R8jJLSAjK6fShOpNSAg0iggjKtJFVGQYjSLDaBQRRmTJh1XY0Q+a8LBQ58cVSpj7eZgrFJcrlLDQEFyuUFyhIbhcIbhCQwgJCSE0NITQ0kcnOYeGhIA7J1aWGotxWhHFxc6LYoopKoKi4mLnp8j5KSzzWERhUTEFhcUUFhY5j0VFFBQUke9+XVBYRH5BEfkFhe5H56egwEk4Ja9LurDyC462RgrLf82vZeFhRz/cS/8twkOJDA9zHiOOfvCXPPd8XZIwKtsvWdN5/8sNZO4/QnzTKMaNCtygLX99gbcE5UV1ZzQwwaUm//1OZDRnRLiLdi2a0K6F9xksjuQVsO9gLvsO5pJ1MJcDh3M5cDiP/YfyOJidx6HsfOcxJ5/sI85PTm4hObkF5OQWAJVfFzNOIo+opEURHhZa2upwWhyeScEzSRzdxnP7MonGY11tTVa8KDmVN+Ydnek/Y19gB2356wugJagqHM+MBiZ41MS/nz9GczaKCKNV8zBaNW/sc0wzP1nLnn05NIttxOjhXejTLb60y+lIntP9lF/amigkv/BoqyO/sIhCd8uksLCYgqIiigqL3S2ZIoqLKdPiKcZpFZVvjRQXF1e4duY0skKcx5KWV0lrLMRppYWGOC21UHeLLTQkhDBXKLv2HkK37iMnt4DoRmH075FAt/ZxFVqAEWFOgnGeH20lOl1eLsLdXVcR4c6y472+V9kXk6GntD6u9zpRlZWhqY1rPr7y1xd4S1DGVCHYPhjKJ8y9+4/w1gKl+UmNqh1PMN3ntyg5lf+tSC39vbKPFLBibTpDT2kdFOc50LeZBOOgLX98gbfxr8ZUIdg+GKpKmNVR8gGckZVDMUc/gBclp9ZgtL6rqd+rvsbj7dpOfR+0ZQnKmCoE2wdDTSXMYPsADrYvAsEWz7iknkSGl53YuCEM2rIuPmOqUJMXg2uiS62mRk8F2wdwsN3WEWzxNNRBW5agjKlCTX0w1NQ1jZpKmMH2ARxst3UEWzzQMAdtWYIy5hhq4oOhpgZb1FTCDLYP4GBrIQRbPA2VJShj/KAmu9RqImEG4wdwsLUQgi2ehsgSlDF+EGxdamAfwCb42Sg+Y/ygoY7CMuZEWAvKGD8Ixi41Y4KdJShj/MS61IypHuviM8YYE5QsQRljjAlKlqCMMcYEJUtQxhhjglJDGCThAkhLSwt0HMYYY8rx+Gx2lV/XEBJUa4Crr7460HEYY4zxrjWw0XNBQ0hQy4HhwC6g8BjbGmOM8S8XTnJaXn5FSHFxccXNjTHGmACzQRLGGGOCkiUoY4wxQckSlDHGmKBkCcoYY0xQsgRljDEmKFmCMsYYE5QsQRljjAlKlqCMMcYEpYYwk8RxE5G+wEtAH2ATMEFVK9ztHEgici7wONAdSAeeVNVpIhIJHATyPDb/VlXPC0CYFYjIBGAakOuxeCLwFjAVuBRn5o9/qepk/0dYlohcjROvpyhgITCaIDvXIjIEmKuqLdyvI6jivIrIZcBjOHf0fwWMV9X0AMfcAngOGAmEAPOB21Q1y71+JnAZUODxNn1UdVMAY67y/12QnudD5TYJAyKBtqq6M5Dn2RKUF+7/0B8CzwJnAWOAz0Sko6oeCGhwbiLSHngXuA4n1oHAAhHZAuwFMlW1VcACrNoA4GlVvctzoYhMBgToCsQCn4rIDlWdGYAYS6nqm8CbJa9FpD/wGfA34FSC5FyLSAhwA/BUuVWT8HJeRaQX8AqQBKwAngBmAyMCHPN0YD/QGQgH3gBeAK5yrx8AXKyqn/ojTk9VxOz1byFYz7OqxnhsEwZ8CSxS1Z3uxQE7z9bF510iEK6qz6pqvqrOBn4GLg9sWGV0Amap6vuqWuRu3S0CzsBJVqsCGNuxeIvvOuBRVc1S1S04/5lu8mdgxyIi4TjJ6iFVXU1wnetJwB+AR8otr+q8XgN8rKpLVPUIcDdwhoh0D1TMIhIKFAGTVPWwqu4D/gOc6V4fBZxM4M67t/Nc1d9C0J3nStyJ82XgQQj8ebYE5V0vYG25Zb/gfEMKCqq6WFVvLnktIs1wJsZdifOtp4WIrBGR3SIyR0TaBipWTyLiwuk2vVZEdorIBhG5S0TicLo+Ujw2D6pz7jYRyAFedL8OpnP9kqoOxPmGDoCINKXq89rLc52qZgOp+O+8V4jZ/YXrYlXd4LHdxTh/2wD9cLqc/iMiGSLyg4hc6Kd4oZKY3ar6Wwi68+xJRNoA9wA3q2qRe3FAz7MlKO9igOxyy7KB6ADEckwiEgt8BHyH0913GPgGp/9ecD5Q3w9YgGUl4PwneR2n++ZSnG92f3Sv9zzvQXXO3V2/f8NpPZXMtBw059qjW8ZTSReOt/Ma0L91LzGXISJ34CSoO92LmgCLcVoFbYBHgbfd141rXRUxV/W3EOzn+S/Ap6rq2VoK6Hm2a1DeHca5CO4pGih/QTHgRKQHTlJKAa52f/u5vdw2twMZItJeVVMDEGYpVU0DzvZYtEpEpuD0zUPZ8x5s5/wCnK6neSULVDVoz7XbYfejt/MatH/r7u7UKTgDUUao6i8AqvoZzjXAEu+KyPXARcBqvwfqVtXfAsF9nl043cBlCucF+jxbC8q7FJxvQJ5Opmw3ScCJyFk4raYPgEvdfduIyMMi0tNj0wj34xE/h1iBiPQWkUnlFkfgxJZG2fMebOf8N8DbHl0gQX2uAdyj3qo6r2X+1kUkGuhAgM+7iDQBPgcGA0M8v9mLyGgRua7cLiV/QwFzjL+FoDzPbsPcjws9Fwb6PFsLyrsvgRAR+QvO8NwxONdNgqWbDBHpCswF7lXVKeVW9wEGiUjJiKfngHmqmuHPGL3YB/xVRLbjjGrqD/wJuBVnIMqDIrIGp0vkDpzYg8VQ4P5yy4L5XJd4A+/ndRawREQSgaXAZGClqq4LRKAeZuN8iR7uvl7jyQU8JyJrgWScwUvDgBv9G2IFXv8WRCRYzzM4f9fLPL94uQX0PFsLygtVzcPpchoDZAL34gy1DKYPnYk4fcSTReSQx88TOMNJs4ANwBac+zKuDVikHlR1B04XwU3AAZyh8v9Q1XeAB4CfcBLVcve6lwIUamU6AeX78oP2XHvwel5V9Udggvv1HqA3MDYwYTpEpA8wChgCpHv8bW8HUNUPcP5PvoXzN/RX4EJV3RaomN28/i0E43n20ImKf9cBP89WUdcYY0xQshaUMcaYoGQJyhhjTFCyBGWMMSYoWYIyxhgTlCxBGWOMCUqWoIwxxgQlu1HXGB+JyAyc6WC8mYQzm/yXQBNV9csUNu5par4Bxh3PTZ8iUgyMVtW5Pmw7FViuqq9XP1JjqsdaUMb47jacWcFb45RjAedG0pJlTwHfup8frmT/2vInYPUJzEjQGmdKIV88DDwsIs2P81jG+Mxu1DXmOIjIKcCPQGd3faVAxdEI2IYzkepPfjrmq8A2VX3IH8czDZd18RlTg9zzrJV28bm7z67EKU4nOGVGrsEp2XEtzvQxd6vqG+79mwBP45QgKQb+h1Pm3FuphCuAfSXJSUQ6AZtxppL6F9AO+AKnnMlTOLOC7wJucc9UXaaLT0QW4ZQi7wuch1Ov6ElVne5xzPeAV0TkUVXNP+6TZcwxWBefMbXvceDPOBNydgB+wElMg3E+7KeJSEnNppdxEtn5OCVJioEF7lLclfk1UFkp7n/glEY/F6c7cg1O9+NAnKJ/0yvZp8SdOF1+/XGS279FxLOE+RdAc/d7GVNrLEEZU/teUNUv3eUi5uLU/7lHVRWnlRMFdBaRLjgtoqtUdbm7VXQtzkSeF3h570E4E8CW96j7PZbgFJxLUdXn3fWUXgDau1trlVmkqi+447sbp6elT8lKd0mXTe5jG1NrLEEZU/s8y5ZnA1s8qvGW1NWJxCkJDqAls3cDe4HGVKxNVqIlzszYxzrmJo/XnsesTOlgC1U94H4aXm6bvUALL/sbUyPsGpQxta/8dZryNXdKhLm37Y/Ttecp08s+RUDICRyzMnmVLCt/DBdQWI33NKbarAVlTPBYi9NSaayqG1R1A86AhieBHl72SQMS/BSfp3j3sY2pNZagjAkS7ms+HwEzRWS4iJwMzMQZXPGLl92ScUbc+Y2IxAIdcQofGlNrLEEZE1yuwxmK/gFOAogFzlXVfV62n4cz2s+fzsRpPa3083FNA2M36hpTh4lINE5p8QtU9Qc/HfMtnFGB//DH8UzDZS0oY+owVc3GuUY10R/HE5HWOC22F/xxPNOwWYIypu57BugjIt6Gotek+4D7VNXbqEJjaox18RljjAlK1oIyxhgTlCxBGWOMCUqWoIwxxgQlS1DGGGOCkiUoY4wxQen/ASvHNQJHt4xeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(results.G, label='simulation')\n",
    "plot(data.glucose, 'bo', label='glucose data')\n",
    "\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration (mg/dL)')\n",
    "\n",
    "savefig('figs/chap08-fig04.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting parameters\n",
    "\n",
    "Based on the parameters of the model, we can estimate glucose effectiveness and insulin sensitivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def indices(params):\n",
    "    \"\"\"Compute glucose effectiveness and insulin sensitivity.\n",
    "    \n",
    "    params: sequence of G0, k1, k2, k3\n",
    "    data: DataFrame with `glucose` and `insulin`\n",
    "    \n",
    "    returns: State object containing S_G and S_I\n",
    "    \"\"\"\n",
    "    G0, k1, k2, k3 = params\n",
    "    return State(S_G=k1, S_I=k3/k2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>S_G</th>\n",
       "      <td>0.0300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>S_I</th>\n",
       "      <td>0.0005</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "S_G    0.0300\n",
       "S_I    0.0005\n",
       "dtype: float64"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indices(best_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under the hood\n",
    "\n",
    "Here's the source code for `run_ode_solver` and `leastsq`, if you'd like to know how they work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def run_ralston(system, slope_func, **options):\n",
      "    \"\"\"Computes a numerical solution to a differential equation.\n",
      "\n",
      "    `system` must contain `init` with initial conditions,\n",
      "     and `t_end` with the end time.\n",
      "\n",
      "     `system` may contain `t_0` to override the default, 0\n",
      "\n",
      "    It can contain any other parameters required by the slope function.\n",
      "\n",
      "    `options` can be ...\n",
      "\n",
      "    system: System object\n",
      "    slope_func: function that computes slopes\n",
      "\n",
      "    returns: TimeFrame\n",
      "    \"\"\"\n",
      "    # the default message if nothing changes\n",
      "    msg = \"The solver successfully reached the end of the integration interval.\"\n",
      "\n",
      "    # get parameters from system\n",
      "    init, t_0, t_end, dt = check_system(system, slope_func)\n",
      "\n",
      "    # make the TimeFrame\n",
      "    frame = TimeFrame(columns=init.index)\n",
      "    frame.row[t_0] = init\n",
      "    ts = linrange(t_0, t_end, dt) * get_units(t_end)\n",
      "\n",
      "    event_func = options.get('events', None)\n",
      "    z1 = np.nan\n",
      "\n",
      "    def project(y1, t1, slopes, dt):\n",
      "        t2 = t1 + dt\n",
      "        y2 = [y + slope * dt for y, slope in zip(y1, slopes)]\n",
      "        return y2, t2\n",
      "\n",
      "    # run the solver\n",
      "    for t1 in ts:\n",
      "        y1 = frame.row[t1]\n",
      "\n",
      "        # evaluate the slopes at the start of the time step\n",
      "        slopes1 = slope_func(y1, t1, system)\n",
      "\n",
      "        # evaluate the slopes at the two-thirds point\n",
      "        y_mid, t_mid = project(y1, t1, slopes1, 2 * dt / 3)\n",
      "        slopes2 = slope_func(y_mid, t_mid, system)\n",
      "\n",
      "        # compute the weighted sum of the slopes\n",
      "        slopes = [(k1 + 3 * k2) / 4 for k1, k2 in zip(slopes1, slopes2)]\n",
      "\n",
      "        # compute the next time stamp\n",
      "        y2, t2 = project(y1, t1, slopes, dt)\n",
      "\n",
      "        # check for a terminating event\n",
      "        if event_func:\n",
      "            z2 = event_func(y2, t2, system)\n",
      "            if z1 * z2 < 0:\n",
      "                scale = magnitude(z1 / (z1 - z2))\n",
      "                y2, t2 = project(y1, t1, slopes, scale * dt)\n",
      "                frame.row[t2] = y2\n",
      "                msg = \"A termination event occurred.\"\n",
      "                break\n",
      "            else:\n",
      "                z1 = z2\n",
      "\n",
      "        # store the results\n",
      "        frame.row[t2] = y2\n",
      "\n",
      "    details = ModSimSeries(dict(success=True, message=msg))\n",
      "    return frame, details\n",
      "\n"
     ]
    }
   ],
   "source": [
    "source_code(run_ode_solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def leastsq(error_func, x0, *args, **options):\n",
      "    \"\"\"Find the parameters that yield the best fit for the data.\n",
      "\n",
      "    `x0` can be a sequence, array, Series, or Params\n",
      "\n",
      "    Positional arguments are passed along to `error_func`.\n",
      "\n",
      "    Keyword arguments are passed to `scipy.optimize.leastsq`\n",
      "\n",
      "    error_func: function that computes a sequence of errors\n",
      "    x0: initial guess for the best parameters\n",
      "    args: passed to error_func\n",
      "    options: passed to leastsq\n",
      "\n",
      "    :returns: Params object with best_params and ModSimSeries with details\n",
      "    \"\"\"\n",
      "    # override `full_output` so we get a message if something goes wrong\n",
      "    options['full_output'] = True\n",
      "\n",
      "    # run leastsq\n",
      "    t = scipy.optimize.leastsq(error_func, x0=x0, args=args, **options)\n",
      "    best_params, cov_x, infodict, mesg, ier = t\n",
      "\n",
      "    # pack the results into a ModSimSeries object\n",
      "    details = ModSimSeries(infodict)\n",
      "    details.set(cov_x=cov_x, mesg=mesg, ier=ier)\n",
      "\n",
      "    # if we got a Params object, we should return a Params object\n",
      "    if isinstance(x0, Params):\n",
      "        best_params = Params(Series(best_params, x0.index))\n",
      "\n",
      "    # return the best parameters and details\n",
      "    return best_params, details\n",
      "\n"
     ]
    }
   ],
   "source": [
    "source_code(leastsq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:** Since we don't expect the first few points to agree, it's probably better not to make them part of the optimization process.  We can ignore them by leaving them out of the `Series` returned by `error_func`.  Modify the last line of `error_func` to return `errors.loc[8:]`, which includes only the elements of the `Series` from `t=8` and up.\n",
    "\n",
    "Does that improve the quality of the fit?  Does it change the best parameters by much?\n",
    "\n",
    "Note: You can read more about this use of `loc` [in the Pandas documentation](https://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-integer)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** How sensitive are the results to the starting guess for the parameters?  If you try different values for the starting guess, do we get the same values for the best parameters?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Related reading:** You might be interested in this article about [people making a DIY artificial pancreas](https://www.bloomberg.com/news/features/2018-08-08/the-250-biohack-that-s-revolutionizing-life-with-diabetes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
